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I.  INTRODUCTION 


This  report  summarizes  three  distinct  studies  performed  during 
the  past  fiscal  year.  The  first  study  was  the  simulation  of  short 
period  regional  phases  for  the  purpose  of  determining  what  physical 
characteri sties  of  such  phases  might  be  used  as  a  seismic 
discriminant.  The  second  study  was  the  expansion  of  the 
capabilities  of  the  three-dimensional  finite  difference  code  TRES 
used  for  three-dimensional  earthquake  simulation.  The  last  study 
was  to  provide  an  algorithm  for  simulation  of  the  near-field  ground 
motion  from  the  1971  San  Fernando,  California,  earthquake  which 
could  be  used  by  the  Mission  Research  Corporation  (MRC)  for  the 

simulation  of  the  acoustic  radiation  generated  by  a  moderate  sized 
event. 

The  regional  phase  synthesis  effort  is  described  in  Sections 
II  through  V.  The  results  reported  here  are  in  addition  to  those 
reported  previously  in  Bache  et  al . ,  (1980),  and  the  results  of  that 
study  are  not  repeated  here.  In  Section  II,  comparisons  are  made  be¬ 
tween  Rayleigh  wave  multimode  synthetic  seismograms  and  those  gener¬ 
ated  by  algorithms  producing  the  complete  layered  media  response.  The 
comparisons  are  quite  favorable  suggesting  that  modal  solutions  are 
adequate  for  describing  the  theoretical  characteristics  of  the  phase 
Lg  at  considerable  savings  compared  to  more  robust  methods. 

In  Section  III,  theoretical  Lg  codas  are  synthesized  using 

earth  models  with  frequency  independent  intrinsic  attenuation.  In 

Bache  et  al . ,  (1980),  it  was  suggested  that  frequency  independent 

models  of  attenuation  were  not  adequate  for  describing  many  of  the 
characteristics  observed  in  Lg.  The  proposed  frequency  dependent 
models  were  found  to  be  no  better  for  describing  observed  behavior 
of  the  short  period  motion.  We  hypothisize  that  many  of  the  common 
features  in  the  data  are  closely  related  to  effects  of  scattering 
due  to  lateral  inhomogenieties  present  along  the  propagation  paths, 
and  that  these  effects  cannot  be  approximated  by  simple  models  of 
intrinsic  attenuation. 
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In  Section  IV,  synthetic  short-period  Lg  codas  are  compared 
for  many  different  crustal  models.  Using  a  published  oceanic 
crustal  model,  synthetic  seismograms  do  not  contain  anything 
resembling  the  phase  Lg.  It  appears  that  the  substantial  velocity 
gradients  prese  t  in  the  oceanic  crustal  model  disperse  the  energy 
uniformly  over  a  very  wide  group  velocity  interval,  and  nothing 
distinctive  is  seen  in  the  window  normally  associated  with  Lg. 
Using  suggested  crustal  models  for  the  Tibetan  Plateau,  synthetic 
seismograms  do  suggest  an  Lg  type  phase.  The  lack  of  an  efficient 
wave-guide  is  not  a  reasonable  explanation  as  to  why  Lg  is  often  not 
seen  in  that  tectonic  province.  For  all  of  the  crustal  models 
considered  here,  it  is  found  that  the  Lg  onset  velocity  is  a  few 
percent  less  than  the  shear  velocity  at  the  base  of  the  crust.  It 
appears  that  the  Lg  coda  begins  with  the  arrival  of  the  shear  wave 
critically  reflected  off  the  Moho. 

In  Section  V,  the  depth  dependence  of  Lg  amplitudes  on  the 
depth  of  dislocation  sources  is  examined.  For  vertical  strike-slip 
and  oblique  thrust  source  mechanisms,  the  amplitudes  of  theoretical 
Lg  is  nearly  independent  of  source  depth  confined  to  the  crust  for 
fixed  mb*  For  vertical  dip-slip  sources,  a  distinct  depth 

dependence  is  seen. 

Two  tasks  under  this  contract,  allocated  approximately  ten 
man-weeks  each,  were  directed  toward  development  of  an  improved 
three-dimensional  seismic  source  modeling  capability,  and  they  are 
described  in  Section  VI.  The  first  of  these  tasks  was  to  provide 
support  to  Technology  Development  Corporation  (TDC)  in  order  to 
modify  the  ILL  I AC  IV  version  of  the  TRES  finite  difference  code. 

The  code  modifications  included  the  following:  1)  a  capability  to 
model  multiple  material  types,  2)  incorporation  of  a  shear-failure 
model  for  earthquake  simulation,  and  3)  a  scheme  for  modeling 
outcropping  thrust  faults. 

The  second  of  these  tasks  was  to  compare  the  performance  of 
the  TRES  code  on  three  computers:  the  ILL  I AC  IV,  S^’s  UNIVAC 

1100/81,  and  the  CRAY  I  computer  at  the  National  Center  for 
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Atmospheric  Research  (NCAR).  This  task  required  conversion  of  the 
UNIVAC  version  of  the  TRES  to  the  NCAR  CRAY,  followed  by  partial 
vectorization  of  the  code  to  exploit  the  CRAY'S  vector  processing 
capabilities.  We  conclude  that  the  TRES  code  now  executes  large 
jobs  with  speed  comparable  to  that  of  the  ILLIAC,  and  we  estimate 


that  full  optimization 

will 

lead  to  further 

improvement 

of 

approximately  a  factor  of 

four 

in  speed. 

Section  VII  describes 

an  algorithm  for 

simulating 

the 

near-field  ground  acceleration  expected  in  the  1971  San  Fernando, 
California,  earthquake.  This  ground  motion  is  to  be  used  as  a 
driver  for  simulating  the  acoustic  radiation  for  that  event.  The 
procedures  established  for  the  simulation  of  the  atmospheric 
disturbances  require  a  relatively  simple  and  efficient  algorithm  for 
determining  the  ground  response  which  is  reasonably  accurate  in  the 
intermediate  frequency  band  (.5  to  3  Hz).  The  source  model  employed 
was  guided  by  previous  studies  of  the  near-field  and  teleseismic 
data.  It  approximates  the  rupture  process  using  three  distinct 
sources  of  radiation.  The  ground  motion  is  computed  using  an 
intuitive  approximation  for  the  radiation  expected  from  a  localized 
propagating  shear  crack.  Comparisons  are  shown  between  the  ground 
velocity  computed  using  the  algorithm  and  the  ground  velocity 
observed  at  Pacoima  Oam. 

Sections  II  through  V  were  prepared  by  T.  C.  Bache  and  H.  J. 
Swanger;  Section  VI  was  prepared  by  S.  M.  Day  and  B.  Shkoller,  and 
Section  VII  was  prepared  by  H.  J.  Swanger. 
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II  COMPARISON  OF  MODAL  AND  TOTAL  SOLUTIONS 


2.1  INTRODUCTION 

The  normal  modes  provide  a  partial  solution  for  the  elastic 
waves  in  a  plane-layered  earth  model.  Our  assumption  is  that  this 
partial  solution  is  a  very  good  approximation  to  the  total  solution 
for  many  problems  in  regional  seismology,  particularly  for  the 
synthesis  of  Lg  and  late  arriving  seismic  energy.  Testing  this 
assumption  requires  comparison  with  an  "exact"  solution. 

Detailed  studies  of  the  performance  of  modal  approximations  of 
SH-motion  have  been  given  by  Herrmann  (1977)  and  Swanger  and  Boore 
(1978).  These  studies  have  demonstrated  that  Love  wave  modes  make 
up  nearly  all  SH  wave  energy  at  distances  beyond  a  few  crustal 
tnickneseses.  The  performance  of  Rayleigh  modes  in  P-SV  problems 
has  not  been  thoroughly  investigated. 

In  our  last  report  (Bache  et  al . ,  1980)  we  described  a 
comparison  of  modal  seismograms  with  seismograms  computed  with  the 
PROSE  program  (Apsel,  1979),  a  direct  wave-number  integration 
program  that  computes  the  total  solution.  This  comparison  resulted 
in  some  ambiguous  results  which  turned  out  to  be  primarily  due  to 
differences  in  the  way  the  filter  representing  the  seismometer 
response  was  specified.  Subsequent  to  the  PROSE  comparisons, 
Henry  Swanger  and  Boris  Shkoller  developed  a  new  wave-number 
integration  program  which  is  similar  in  many  respects  to  PROSE. 
This  program  was  used  to  compute  total  seismograms  for  comparison 
with  the  modal  solution.  The  comparison  confirms  that  the  P-SV 
modes  give  an  excellent  approximation  to  the  total  solution  for 
arrivals  with  group  velocities  less  than  the  shear  velocity  of  the 
underlying  halfspace. 

2.2  COMPARISON  WITH  PROSE 

The  PROSE  calculations  were  done  by  Dr.  John  Orcutt  of  the 
Scripps  Institute  of  Oceanography.  The  structure  is  listed  in  Table 
1  and  the  phase  and  group  velocity  dispersion  for  all  the  modes  (66) 
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TABLE  1 


CRUSTAL  MODEL  FOR  COMPARING  COMPLETE  AND  MODAL  SEISMOGRAMS 


Oepth 

(km) 

Thickness 

(km) 

a 

(km/sec) 

8 

(km/sec) 

(gm/cm^) 

0 

1.5 

1.5 

3.7 

2.16 

2.10 

35 

8.0 

6.5 

6.1 

3.3 

2.85 

250 

34.0 

26.0 

6.6 

3.59 

3.05 

1800 

00 

00 

8.1 

4.52 

3.35 

2000 

5 


between  0  and  5  Hz  is  plotted  in  Figure  1.  The  source  was  a  strike- 

slip  double-couple  at  a  depth  of  1  km.  The  azimuth  was  45°  from  the 

strike  and  the  source  time  function  was  a  step  with  a  moment  of 
22 

10  dyne-cm.  The  WWSSN  short  period  instrument  response  was 
included  and  the  spectrum  was  tapered  to  zero  between  0.5  and  1.0 
Hertz  with  a  cosine-squared  filter.  Thus,  the  seismograms  include 
no  energy  above  1  Hertz  (only  14  modes  contribute  to  the  modal 
solution) . 

The  comparison  of  the  modal  and  PROSE  seismograms  from  Bache, 
et  al_,  (1980)  is  shown  in  Figure  2.  The  waveform  agreement  is  quite 
good  and  indicates  that  most  of  the  important  energy  is  included  in 
the  modal  solution.  The  amplitude  comparison  is  not  as  good.  This 
turned  out  to  be  primarily  due  to  a  different  specification  of  the 
WWSSN  short  period  seismometer  response.  More  complete  comparisons 
which  are  done  in  an  entirely  consistent  way  are  described  in  the 
next  section. 


2.3  COMPARISON  WITH  S3,S  DIRECT  WAVE  NUMBER  INTEGRATION  PROGRAM 

A  new  direct  wave-number  integration  program  was  developed  by 

3 

Henry  Swanger  and  Boris  Shkoller  of  S  to  compute  the  complete 
solution.  This  program  does  essentially  the  same  calculation  as 
PROSE,  but  the  formulation  and  numerical  procedures  are  different  in 
detai 1 . 

The  first  comparison  with  the  results  of  this  new  direct 
wave-number  integration  program  is  shown  in  Figure  3.  The 
seismograms  were  computed  the  same  as  the  R  =  50  km  case  in  Figure 
2,  except  that  the  structure  (Table  1)  was  modified  to  include 
infinite  Q. 

The  complete  solution  was  done  with  a  Nyquist  frequency  of  1.0 
Hz,  hence  the  time  step  was  0.5  seconds  and  the  seismograms  have  a 
jagged  appearance.  The  time  series  for  the  modal  calculations  was 
sine  interpolated  to  a  time  step  of  0.1  seconds  by  filling  the 
spectrum  with  zeros  between  1  Hertz  and  a  5  Hertz  Nyquist  frequency. 
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3.5  2.0 


12.5  25.0  37.5  50.0  52.5  75.0  37.5  100.0  112.5 


TIME  (s) 


Figure  2.  Synthetic  seismograms  computed  with  modal  superposition 
are  compared  with  those  computed  with  PROSE.  Zero  on 
the  time  scale  of  the  plots  is  reduced  by  R/8.1.  The 
maximum  peak-to-peak  amplitude  in  microns  is  shown  at 
the  left  and  several  interesting  group  velocities  are 
marked  in  each  set. 
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Several  more  comparisons  are  shown  in  Figures  4  to  7.  First, 
in  Figure  4  we  repeat  the  comparison  of  Figure  3,  this  time 

including  the  Q.  The  same  comparison  is  shown  in  Figure  5  without 
the  WWSSN  short  period  instrument  response.  In  Figure  6  the 

seismograms  (with  Q  and  the  seismometer  included)  are  shown  for  a 
source  depth  of  3.0  km. 

The  most  complete  comparison  is  shown  in  Figure  7.  The  source 
and  azimuth  are  the  same  as  for  all  comparisons  in  this  section. 
The  depth  is  3.0  km  and  the  elastic  (infinite  Q)  version  of  the 

structure  in  Table  1  was  used.  This  time  the  calculations  were  done 
with  a  Nyquist  frequency  of  5.0  Hz.  The  phase  and  group  velocity 
dispersion  for  the  modes  were  shown  in  Figure  1.  There  are  66  modes 
from  0  to  5  Hertz.  The  seismograms  in  Figure  7  include  the  WWSSN 
short  period  instrument  response  and  the  spectra  for  both  cases  were 
tapered  between  4.5  and  5.0  Hertz  with  a  cosine-squared  filter. 

The  comparisons  in  Figure  3  to  7  demonstrate  that  the  modes 

give  a  remarkably  accurate  approximation  for  the  complete  ground 
motion  in  the  group  velocity  window  <  4.5  km/sec.  The  comparison  in 
Figure  7  is  a  particularly  good  demonstration  of  this  point.  The 
complete  and  modal  solutions  have  almost  the  same  phasing;  that  is, 
zero  crossings  and  breaks  in  the  waveform  occur  at  nearly  the  same 

time  in  both.  Differences  between  the  two  are  mainly  in  the 

amplitudes  of  the  peaks,  but  even  these  differences  are,  at  most,  30 
to  40  percent. 

Our  comparison  suggests  that  Rayleigh  modes  are  adequate  for 
describing  the  characteristics  of  Lg.  In  this  study,  we  will 
restrict  our  attention  to  modal  synthesis  of  Lg  even  though 
complete  solutions  are  obtainable.  Modal  solutions  have  two  major 
advantages  over  complete  solutions.  First,  modal  solutions  are 

constructed  from  intermediate  results  which  allow  for  better 
interpretation  of  the  final  seismograms.  As  demonstrated  in  our 

last  report  (Bache  et  aj^,  1980),  the  dispersion  parameters  and 

modal  energies  provide  a  great  deal  of  information  lost  within  the 
"black-box"  of  complete  seismogram  programs.  Second,  modal 
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7.39 


Figure 


Figure  5.  The  comparison  in 


gure 


Figure  7.  Seismoorams  like  those  in  Fiqure  6  are  compared.  The  differences  are  that  the 
infinite  and  that  the  calculations  are  carried  out  to  5  Hz. 


solutions  have  significant  cost  advantages  over  complete  solutions. 
Our  experiences  suggest  that  the  cost  of  obtaining  all  the  needed 
dispersion  parameters  for  a  given  earth  structure  is  three  or  four 
times  less  than  that  needed  to  compute  a  single  complete  regional 
seismogram.  Once  the  dispersion  functions  are  obtained,  a  modal 
seismogram  costs  about  one-hundredth  that  of  a  complete  seismogram. 
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III.  SYNTHESIS  OF  Lq  WITH  FREQUENCY-DEPENDENT  Q 


3.1  INTRODUCTION 

In  our  semi-annual  report  (Bache  et  al_. ,  1980)  we  focused  our 
attention  on  the  synthesis  of  Lg  using  plane-layered,  laterally 
homoqeneous  earth  models  with  frequency-independent  0.  While  many 
of  the  important  qualitative  features  of  observed  Lg  can  be 
reproduced  with  such  models,  we  found  there  were  three  basic 
deficiencies  that  appeared  to  be  impossible  to  remove,  as  long  as  we 
are  confined  to  this  restricted  class  of  earth  models.  These  are: 

1.  The  synthetics  do  not  include  the  high 

frequency  enerqy  arriving  with  apparent 

velocities  slower  than  2.9  km/sec  that  is  a 
prominent  feature  of  the  observations. 

2.  The  synthetic  spectra  fall  off  more  rapidlv  at 
high  frequencies  than  the  observed  spectra  from 
the  SALMON  nuclear  explosion. 

3.  The  amplitude  attenuation  with  range  is  much 
faster  than  is  observed  in  the  eastern  United 
States  for  earthquakes  or  explosions. 

To  what  can  we  attribute  these  deficiencies?  The  failure  to 
include  multi pathinq  and  scattering  by  lateral  heterogeneities  is 
likely  to  contribute  to  the  first  of  them.  But  the  other  two  (and, 
perhaps,  some  of  the  first)  are  likely  to  be  caused  by  the 
assumption  of  frequency-independent  Q. 

In  this  section  we  study  the  effects  of  allowing  the  Q  to  be 
frequency-dependent.  In  evaluating  the  models  we  are  mostly 
concerned  with  their  ability  to  replicate  the  observed  amplitude 
attenuation.  This  is  observed  for  many  events  at  different  depths 
and  locations.  The  spectral  characteristics  are  also  important,  but 
are  more  difficult  to  evaluate  because  of  uncertainties  about  the 
source. 


16 


3.2  SYNTHETIC  SEISMOGRAMS  WITH  THE  MITCHELL  (1980)  VARIABLE  Q  MODEL 

Mitchell  (1980)  points  out  that  frequency-independent  Q  models 
inferred  from  fundamental  mode  Rayleigh  wave  observations  in  eastern 
North  America  are  not  consistent  with  the  observed  attenuation  of  Lg 
in  this  region.  This  is,  of  course,  one  of  the  points  made  in  our 
semi-annual  report.  To  resolve  this,  Mitchell  (1980)  assumed 

Q^fu.z)  =  C(z)  u'C,  (1) 

D 

and  obtained  several  Q”*(u,z)  models  corresponding  to  different 

D 

Observations  of  the  attenuation  of  Lg  between  1  and  10  seconds 
in  eastern  North  America  were  then  used  to  select  the  best  values 
for  l,  . 

Assuming  Q  is  frequency-independent,  Mitchell  (personal 
comnuni cation)  finds  that  the  best  fit  to  the  data  is  provided  by 
either  of  the  models  Sll  or  S12  listed  in  Table  2.  The  velocity 
profile  is  based  on  the  model  fit  to  Rayleigh  waves  by  McEvilly 
(1964),  which  has  a  38  km  crust.  This  model,  which  is  discussed  in 
more  detail  in  Section  4.4,  has  been  modified  by  adding  some  lower 
velocity  material  in  the  top  4.1  kilometers.  As  was  discussed  by 
Bache,  et  aj_.  (1980),  this  seems  to  be  necessary  to  avoid  having  the 
short  period  synthetics  be  dominated  by  a  pulse-like  arrival 
associated  with  the  fundamental  mode  trapped  in  the  near-surface 
1 ayers. 

Synthetic  seismograms  for  models  Sll  and  S12  are  plotted  at 
two  ranges,  500  and  1000  km,  in  Figure  8.  Fifty  modes  were  used  for 
these  calculations.  The  source  is  a  reduced  displacement  potential 
for  the  SALMON  explosion  given  by  Murphy  (1969)  and  is  at  a  depth  of 
0.83  kilometers,  which  is  appropriate  for  SALMON.  The  LRSM  short 
period  seismometer  response  is  included  in  the  synthetics.  The 
synthetics  for  the  two  models  are  very  similar,  with  the  main 
difference  being  in  the  peak  amplitudes. 

The  SALMON  observations  recorded  by  LRSM  seismometers  at 
ranges  near  250  and  1200  km  are  plotted  in  Figure  9  at  the  same  time 
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TABLE  2 


FREQUENCY-DEPENDENT  Q  MODELS  FOR  THE  EASTERN  U.S. 
(Mitchell,  personal  communication) 


Sll  S12 
(5-0.3)  (5-0.2) 


Layer 

Depth 

(km) 

Thickness 

(km) 

a 

(km/sec ) 

8 

(km/sec) 

p 

(km/sec ) 

C 

C 

1 

0.6 

0.6 

3.7 

2.16 

2.1 

245 

245 

2 

1.0 

0.4 

4.55 

2.54 

2.2 

575 

480 

3. 

2.6 

1.5 

4.55 

2.54 

2.2 

610 

505 

4 

4.1 

1.5 

5.60 

3.14 

2.65 

590 

510 

5 

11.0 

6.9 

6.10 

3.50 

2.70 

500 

410 

6 

16.0 

5.0 

6.40 

3.68 

2.90 

630 

490 

7 

20.0 

4.0 

6.40 

3.68 

2.90 

990 

750 

8 

25.0 

5.0 

6.70 

3.67 

2.90 

1400 

1040 

9 

31.0 

6.0 

6.70 

3.67 

2.90 

2280 

1700 

10 

38.0 

7.0 

6.70 

3.67 

2.90 

2780 

2160 

11 

43.0 

5.0 

8.15 

4.67 

3.30 

3235 

2600 

12 

49.0 

6.0 

8.15 

4.67 

3.30 

5075 

4170 

13 

oo 

oo 

8.15 

4.67 

3.30 

6150 

5200 
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500  km 


Figure  8.  Synthetic  seismograms  for  SALMON  are  plotted  for  two  crustal  models  at  two  ranges.  The 

source  is  a  reduced  displacement  potential  from  Murphy  (1969)  at  a  depth  of  0.83  kilometers. 
The  LRSM  short  period  seismometer  response  is  included. 
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Figure  9.  Observed  short  period  recordings  of  SALfN 
and  minimum  amplitudes  are  indicated  aloi 
event  origin  time  was  16:00.0).  The  timi 
indicated  on  each  record. 
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lings  of  SALMON  are  shown  from  eight  LRSM  stations.  The  maximum 
Indicated  along  with  the  time  of  zero  on  the  time  scale  (the 
10).  The  times  for  several  apparent  group  velocities  are 
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scale  as  the  synthetics  In  Figure  8.  Comparing  the  two,  it  appears 
that  the  synthetics  have  too  much  high  frequency  energy,  especially 
at  the  larger  range.  Also,  there  is  no  indication  in  the  data  of 
the  large  amplitudes  that  arrive  with  an  apparent  velocity  of  about 
3.05  km/sec  on  the  R  =  500  km  synthetics.  This  arrival  can  also  be 

seen  at  R  =  1000  km,  though  it  is  less  prominent. 

For  evaluating  synthetics  like  those  in  Figure  8,  we  will 

primarily  examine  the  apparent  amplitude  attenuation  of  Lg.  This  is 
a  concept  introduced  by  Nuttli  (1973),  who  assumed  that  Lg 
attenuates  with  distance  as 

sinA-1^2  exp(-rA),  (2) 

where  the  a  is  the  range  and  y  is  the  apparent  attenuation. 

Plotting  Lg  amplitude  versus  range,  Nuttli  (1973)  found  the  y 
providing  the  best  fit  to  the  trend  of  the  data.  This  technique  has 
been  used  by  several  others  to  fit  various  sets  of  Lg  amplitudes  in 
eastern  North  America  (e.g.,  Bollinger,  1973,  1979;  Street,  1976). 
The  consensus  is  that  the  y  for  this  region  is  about  0.07  to  0.11 

deg-1.  The  sparsely  sampled  Lg  amplitudes  for  the  SALMON  data  in 

Figure  9  are  consistent  with  these  values,  though  they  are  probably 
best  fit  with  slightly  larger  y,  about  0.12  deg-1  or  so  (see 

Figure  10). 

The  convention  for  the  Lg  amplitude  on  the  synthetics  and 

SALMON  data  was  to  measure  the  maximum  (3  cycles  or  so)  sustained 
amplitude  within  about  ten  seconds  of  the  arrival  time  associated 
with  a  3.5  km/sec  group  velocity.  The  synthetic  y  is  then  computed 
from  (2).  That  is, 


where  A  is  the  Lg  amplitude  and  subscripts  1  and  2  refer  to  the  two 
ranges. 
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The  Lg  amplitudes  for  the  two  models  are  listed  in  Table  3. 
Also  listed  are  Lg  amplitudes  from  the  SALMON  observations  in  Figure 
9.  Using  (3),  the  y$yn  for  both  models  is  0.18  deg"1,  which  is 
larger  than  the  observed  y.  Also,  the  Lg  amplitudes  are  very  much 
larger  than  the  observed  values  for  SALMON.  All  these  amplitude 
data  are  plotted  in  Figure  10. 

Do  the  models  give  a  y  more  like  the  data  for  deeper 

earthquake  sources?  To  address  this  question  we  computed  the 

synthetics  shown  in  Figure  11.  In  this  case,  the  source  is  a 

strike-slip  point  double-couple  at  a  depth  of  5  km,  with  a  step 

25 

function  time  history  represented  by  a  moment  of  10  dyne-cm. 
The  synthetics  are  computed  at  an  azimuth  30®  from  the  strike  and 
the  WWSSN  short  period  seismometer  response  is  included.  As  with 
the  explosion  source,  the  two  models  give  nearly  the  same  synthetics 
except  for  an  amplitude  scaling.  Measuring  the  Lg  amplitudes  on 
those  synthetics  and  computing  y  from  (3),  we  find  that  ysyn 
=  0.20  deg'1  for  both  models.  Thus,  we  see  that  the  y$yn  is  not 
very  dependent  on  the  source  type. 

In  his  analysis  leading  to  the  frequency-dependent  Q  models  in 
Table  3,  Mitchell  (1980)  considers  the  excitation  functions  of  the 
first  six  higher  modes  and  concludes  that  the  results  of  his  study 
do  not  depend  critically  on  the  number  of  modes  employed  at  any 
period.  To  check  the  validity  of  the  assumption,  we  repeated  the 
SALMON  synthetics  for  model  S12  from  Figure  8,  this  time  with  seven 
modes  (rather  than  the  fifty  used  before).  The  resulting  synthetic 
seismograms  are  shown  in  Figure  12. 

The  seven  mode  synthetic  seismograms  are  not  very  realistic 
looking  when  filtered  to  emphasize  the  short  periods,  as  they  are  in 
Figure  12.  These  seismograms  are  dominated  by  the  fundamental 
mode.  It  is  yet  higher  modes,  which  are  associated  with  propagation 
at  the  base  of  the  crust,  that  are  most  important  for  the  high 
frequency  energy  near  3.5  km/sec.  This  is  discussed  in  more  detail 
by  Bache,  et  al .  ( 1980) . 
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TABLE  3 


THEORETICAL  Lg  AMPLITUDES  FROM  FIGURE  8  AND 
OBSERVED  SALMOi.  Lg  AMPLITUDES 


Model 

Range 

(km) 

Lg* 

(nm) 

SI1 

500 

6950 

Sll 

1000 

1710 

S12 

500 

5310 

S12 

1000 

1340 

SALMON  OBSERVED 

Station 

JELA 

244 

2746 

EUAL 

245 

2649 

BLWV 

1057 

259 

VOIO 

1251 

279 

BRPA 

1374 

240 

WFMN 

1427 

153 

*  Uncorrected  for  frequency-dependent  seismometer  response. 
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Figure  12.  The  SALMON  synthetics  for  Model  S12  are  shown  when  computed  with  just  seven  modes. 
Note  that  the  time  scale  is  different  on  the  two  plots. 


The  Lg  amplitudes  for  the  seismograms  in  Figure  12  are  1430  nm 
at  500  km  and  240  nm  at  R  =  1000  km.  These  are  much  smaller  than 
the  amplitudes  of  5310  nm  and  1340  nm  from  the  fifty  mode  synthetics 
in  Figure  8  (Table  3).  Using  (3),  the  Ysyn  for  the  seven  mode 
synthetics  is  0.27.  This  is  higher  than  the  y  for  the  50  mode 
synthetics  because  the  lower  modes  are  associated  with  propagation 
at  shallower  depths  where  0  is  lower. 

We  are  forced  to  the  conclusion  that  there  is  no  simple  way  to 
characterize  the  attenuation  of  Lg  using  qualitative  ideas,  as  in 
Mitchell  (1980).  It  is  necessary  to  consider  a  large  number  of 
modes,  which  is  probably  only  practical  by  examining  synthetic 
seismograms.  From  the  synthetic  Lg  seismograms  we  find  that  the  Lg 
amplitude  attenuation  is  much  greater  than  that  observed. 
Therefore,  the  frequency-dependent  model  proposed  by  Mitchell  (1980) 
is  not  consistent  with  observed  Lg  amplitude  attenuation  in  Eastern 
North  America,  at  least  in  the  band  passed  by  the  LRSM  and  WWSSN 
short  period  seismometers. 

3.3  SYNTHETIC  Lg  SEISMOGRAMS  WITH  MODIFIED  VERSION  OF  MITCHELL'S 

(1980)  FREQUENCY-DEPENDENT  Q  MODEL 

3.3.1  Introduction 

The  frequency-dependent  Q  models  proposed  by  Mitchell  (1980) 
fail  to  give  Lg  seismograms  that  match  essentially  the  same  three 
data  features  listed  in  the  introduction  (Section  3.1)  to  this 
section.  That  is,  the  synthetics  do  not  match  the  duration  beyond 
2.9  km/sec,  the  spectral  content  and  the  apparent  amplitude 
attenuation.  The  frequency-dependent  Q  model  has  very  little  effect 
on  the  duration.  This  is  an  indication  that  the  observed  energy 
arriving  later  than  2.9  km/sec  is  probably  more  strongly  associated 
with  multi-pathing  and  scattering  by  lateral  heterogeneities  than 
with  Q. 
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The  frequency-dependent  Q  model  does  allow  much  more  high 
frequency  energy  into  the  Lg  spectrum.  In  fact,  it  changes  the 
spectral  content  too  far  in  this  direction;  the  synthetics  now 
appear  to  have  too  much  high  frequency  energy.  Finally,  we  have 
apparent  amplitude  attenuation  (rsyn),  the  observed  data  feature 
that  we  are  most  concerned  with  reproducing  in  this  section.  The 
frequency-dependent  models  Sll  and  S12  have  too  large  a  YSyn*  At 
the  same  time,  the  absolute  amplitudes  are  much  larger  than  those 
observed  for  SALMON,  assuming  the  SALMON  source  is  reasonably  well 
known.  Thus,  the  model  we  are  seeking  must  have  more  attenuation  to 
reduce  the  absolute  amplitudes  and  the  high  frequency  content  of  the 
synthetic  Lg  phases.  At  the  same  time  it  must  have  a  lower  YSyn> 
which  suggests  higher  Q.  These  requirements  appear  to  be 
contradictory,  and  so  will  be  very  difficult  to  fulfill.  In  this 
section  we  describe  some  experiments  to  see  just  how  difficult. 

3.3.2  Variations  in  the  Q  Near  the  Surface 

The  first  experiment  is  to  change  the  Q  in  the  top  layers. 
The  Q  in  these  layers  is  not  very  well  resolved  by  Mitchell's 
procedure,  so  we  should  have  some  flexibility.  The  first  three 
models  are  called  S13,  S14  and  S15.  The  Q  for  these  models,  and  for 
Sll  and  S12  from  Table  2,  are  listed  in  Table  4.  Only  the  Q  in  top 
4.1  kilometers  is  altered. 

Synthetic  seismograms  like  those  in  Figure  8  are  plotted  for 
models  S13,  S14  and  S15  in  Figure  13.  The  source  is  that  for 
SALMON.  The  waveforms  are  little  different  from  those  in  Figure  8. 
The  Lg  amplitudes  are  summarized  in  Table  5  (the  model  S16  will  be 
discussed  later),  together  with  the  ^Syn  computed  from  (3).  We 
see  that  lowering  the  Q  near  the  surface  has  had  the  effect  of 
increasing  the  effective  attenuation  of  the  entire  Lg  phase.  The 
amplitudes  are  smaller  and  the  Tsyn  is  larger.  We  conclude  that 
no  minor  fixes,  like  those  represented  by  the  models  in  Table  5, 
will  make  these  frequency-dependent  Q  models  compatible  with  the 
three  important  characteristics  of  the  data  that  we  have  been 
discussing. 
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TABLE  4 


MODIFIED  FREQUENCY  DEPENDENT  Q  MODELS 
FOR  THE  EASTERN  UNITED  STATES 


Layer 

Depth 

(km) 

8 

(km/sec) 

S12 

(5=0.2) 

S13 

(5=0.2) 

Sll 

(5=0.3) 

S14 

(5=0.3) 

S15 

(5=0.: 

1 

0.6 

2.16 

245 

50 

245 

50 

25 

2 

1.0 

2.54 

480 

100 

575 

250 

125 

3 

2.6 

2.54 

505 

505 

610 

300 

150 

4 

4.1 

3.14 

510  _ 

510 

590 

300 

150 

S12 

Model 

Sll 

Model 
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TABLE  5 


Lg  AMPLITUDES  FROM  SYNTHETIC  SEISMOGRAMS 


Model 

R=500 

Lg 

(nm) 

R=1000 

Lg 

(nm) 

’fsyn 

Sll 

6950 

1710 

0.18 

S12 

5310 

1340 

0.18 

S13 

3070 

570 

0.25 

S14 

3730 

790 

0.22 

S15 

2460 

343 

0.31 

S16 

3480 

833 

0.19 

f 


3.3.3  A  Modified  Velocity  Model 

The  Q  inversions  described  by  Mitchell  (1980)  were  done  with 
the  McEvilly  (1964)  model  which  best  fit  his  Rayleigh  wave  data.  We 
have  modified  this  model  by  lowering  the  velocities  in  the  top  4.1 
kilometers,  but  this  is  a  region  that  is  poorly  resolved.  McEvilly 
(1964)  also  gives  a  model  which  gives  the  best  simultaneous  fit  to 
both  Rayleigh  and  Love  wave  data.  This  model  differs  from  that  in 
Table  3  (except  for  the  top  4  layers  that  we  have  introduced)  by 
having  the  shear  velocity  at  the  base  of  the  crust  be  3.94  km/sec, 
rather  than  3.67  km/sec.  This  has  the  effect  of  introducing  some 
gradient  at  the  base  of  the  crust,  rather  than  having  a  27  km 
section  of  nearly  constant  shear  velocity. 

A  crustal  model  with  the  frequency-dependent  Q  of  S14  and  b  = 
3.94  km/sec  between  20  and  38  km  depths,  is  called  S16.  Seismograms 
for  this  model  are  plotted  in  Figure  14  and  the  amplitudes  are 
listed  in  Table  5. 

The  seismograms  for  this  model  have  the  energy  more  evenly 
spread  through  the  Lg  time  window,  especially  at  R  =  500  km.  They 
also  have  a  y  more  like  that  of  Sll  and  S12,  while  the  amplitudes 
are  much  smaller. 

3.4  CONCLUSIONS 

In  our  previous  report  (Bache  et  aj_. ,  1980),  it  was  suggested 
that  frequency  independent  Q  models  are  not  sufficient  to  describe 
the  attenuation  character istics  of  the  observed  Lg  from  SALMON.  In 
this  section  we  have  examined  the  characteristics  of  synthetic  Lg  in 
EUS  models  with  frequency  dependent  Q.  Mitchell  (1980)  has 
suggested  intrinsic  attenuation  with  frequency  power  laws.  For 
models  suggested  from  inversion  studies,  synthetic  Lg  were  found  to 
have  the  same  deficiencies  as  those  with  frequency  independent  Q. 
Modifying  the  properties  of  the  near  surface  layering  gave  some 
improvement,  but  the  result  synthetic  seismograms  are  not  more 
realistic  than  those  found  for  frequency  independent  Q  models. 
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The  major  deficiencies  in  the  synthetic  Lg  -  1)  too  little 
late  arriving  energy,  2)  wrong  spectra  above  3  Hz,  and  3)  wrong 
attenuation  rate  of  peak  amplitudes  with  distance  -  suggest  that 
simple  attenuation  mechanisms  in  plane  layered  earth  models  is  not 
realistic  for  the  high  frequencies  and  large  distances  considered 
here.  The  presence  of  late  arriving  energy  and  lack  of  a 
fundamental  mode  in  the  short  period  data  suggests  that  scattering 
is  probably  an  important  contribution  to  the  attenuation  process. 
Frequency  dependent  Q  models  do  not  appear  to  approximate  the 
effects  of  such  mechanisms  very  well.  Several  models  for  the 
scattering  of  high  frequencies  have  been  suggested  (Chernov,  1960; 
Aki  1969;  Dainty  and  Toksoz,  1977)  and  these  models  are  being 
applied  to  the  study  of  near-regional  seismic  codas  (Aki  and  Chouet, 
1975,  for  example).  Such  models  of  scattered  wave  propagation  might 
be  useful  to  regional  studies,  but  at  this  time  it  is  not  clear  just 
now  they  may  be  applied  to  synthetic  seismogram  computations. 
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IV.  SYNTHETIC  lg  SEISMOGRAMS  IN  DIFFERENT  TECTONIC  PROVINCES 


4.1  INTRODUCTION 

An  important  feature  of  Lg  is  that  its  characteristics  are 
observed  to  vary  widely  from  one  tectonic  province  to  another.  The 
fact  that  Lg  does  not  propagate  on  paths  that  cross  the  deep  ocean 
was  noted  in  the  earliest  papers  describing  this  phase  (Press  and 
Ewing,  1952).  Recent  work  has  concluded  that  the  Lg  phase  is  also 
absent  for  paths  crossing  the  Tibetan  Plateau  (Ruzaiken,  et  al . , 
1977),  and  the  Black  Sea  (Antonova,  et  aj_. ,  1978). 

Several  explanations  for  the  relative  efficiency  of 
propagation  of  Lg  in  different  tectonic  provinces  have  been 
proposed.  There  are  basically  three  reasons  why  the  Lg  may  be 
absent.  Tnese  are: 

1.  The  waveguide  for  Lg  is  absent. 

2.  The  phase  is  disrupted  by  passage  across  a 
sharp  lateral  discontinuity. 

3.  The  phase  is  strongly  attenuated. 

In  this  section  we  study  the  relative  excitation  of  Lg  in 

three  different  tectonic  regions  using  a  straightforward  technique. 
Tne  regions  are  the  deep  ocean,  the  Tibetan  Plateau  and  the  stable 
platform  area  of  the  central  and  eastern  United  States.  The 

technique  is  to  compute  Lg  seismograms  in  standard  models  for  each 
region  and  compare  them.  By  standard  models  we  mean  laterally 
homogeneous,  plane-layered  crustal  models  with  frequency-independent 
Q.  We  found  in  Section  III  that  the  introduction  of  frequency- 
dependent  Q  did  not  change  the  basic  character  of  the  synthetic  Lg 

seismograms  in  terms  of  the  time  of  Lg  onset  or  its  duration. 

With  the  straightforward  technique  we  are  using  we  can 
investigate  the  first  and  third  listed  reasons  for  the  failure  of  Lg 
to  propagate  in  certain  regions.  An  analysis  of  the  effect  of 
lateral  boundaries  would  require  techniques  like  Gregerson  and  Alsop 
(1974)  developed  for  propagating  the  Lg  phase  across  such 
boundaries. 


Synthetic  seismograms  for  the  oceanic,  Tibetan  Plateau,  and 
central  and  eastern  United  States  models  are  presented  in  Sections 
4.2,  4.3  and  4.4,  respectively.  These  seismograms  are  compared  in 
Section  4.5  where  our  conclusions  are  summarized. 

4.2  Lg  IN  AN  OCEANIC  CRUST 

An  often  mentioned  characteristic  of  Lg  is  that  it  is  not 
observed  when  the  travel  path  includes  a  section  that  crosses  a  deep 
ocean  basin  (Press  and  Ewing,  1952).  Several  reasons  have  been 
proposed  to  explain  the  disappearance  of  Lg  on  such  paths.  Some 
have  suggested  that  the  crustal  waveguide  in  which  the  Lg  energy  is 
trapped  is  disrupted  by  the  ocean-continent  margin  (Ewing,  Jardetzky 
and  Press,  1957).  Knopoff,  et  al_.  (1979)  suggest  that  Lg  is 
scattered  by  strong  lateral  variations  in  the  relatively  thick  low 
velocity  sediments  at  the  ocean  bottom  or  is  rapidly  attenuated  by 
the  low  Q  which  is  characteristic  of  the  oceanic  crust. 

In  this  section  we  look  at  the  propagation  of  Lg  in  oceanic 
structures  in  a  straightforward  way.  We  simply  compute  multimode 
seismograms  for  a  good  model  for  the  oceanic  crust  and  see  if  Lg  is 
present.  The  model  is  plane-layered  and  laterally  homogeneous  and 
has  a  frequency-independent  Q.  Models  of  this  kind  for  the 
continental  crust  lead  to  seismograms  with  an  obvious  Lg  phase,  as 
is  demonstrated  elsewhere  in  this  report. 

The  oceanic  crustal  model  is  listed  in  Table  6.  This  model  is 
a  smoothed  version  of  the  model  FF2  obtained  by  Spudich  and  Orcutt 
(1980)  by  fitting  synthetic  seismograms  to  the  data  from  a 
refraction  line  east  of  Guadaloupe  Island  off  the  coast  of  Baja 
California.  We  also  modified  the  model  by  adding  a  higher  velocity 
layer  for  the  half space  below  a  depth  of  19  km. 

There  are  46  modes  for  this  model  between  1  and  5  Hertz.  They 
are  plotted  in  Figure  15,  The  first  thing  to  be  noticed  is  that 
there  is  no  clustering  of  group  velocity  stationary  phases  near  3.5 
km/sec,  as  is  prominently  seen  on  similar  plots  for  continental 
crustal  structures.  The  only  prominent  group  velocity  plateaus  are 
near  1.0  km/sec  and  0.25  km/sec. 
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TABLE  6 


OCEANIC  CRUSTAL  MODEL  FF2  (SPUDICH  AND  ORCUTT,  1980) 


Depth 

T.-.ickness 

a 

s 

p 

Qa 

Qb 

(km) 

(km) 

(km/sec) 

(km/sec) 

(km/sec) 

3.40 

3.40 

1.463 

0 

1.0 

00 

0 

3.57 

0.17 

1.58 

0.25 

1.5 

200 

100 

4.00 

0.43 

4.8 

2.45 

2.07 

450 

100 

4.54 

0.54 

6.0 

3.25 

2.53 

450 

100 

5.12 

0.58 

6.3 

3.76 

2.64 

450 

225 

8.64 

3.52 

6.95 

3.74 

2.89 

450 

225 

9.33 

0.69 

7.38 

4.3 

3.05 

450 

225 

19.0 

9.67 

7.80 

4.5 

3.20 

450 

225 

8.20 

4.62 

3.35 

450 

225 
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Synthetic  seismograms  were  computed  in  the  oceanic  model  for  a 

strike-slip  double-couple  source.  The  source  time  history  was  a 

25 

step  with  a  moment  of  10  dyne-cm.  The  synthetics  were  computed 
for  depths  of  8.4  km  (5  km  below  the  ocean  bottom),  10.9  km  and  12.4 
km.  The  range  was  200  km.  We  also  computed  a  synthetic  at  R  =  100 
km  for  the  shallowest  source.  These  synthetic  seismograms  are  shown 
in  Figure  16. 

The  seismograms  in  Figure  16  have  the  energy  distributed 
rather  evenly  over  the  entire  group  velocity  interval  from  4.62 
km/sec  (the  cutoff  velocity)  to  less  than  2.0  km/sec.  There  is  no 
sharp  arrival  of  energy  near  3.5  km/ sec  which  characterizes  Lg.  We 
will  have  more  to  say  about  these  seismograms  in  Section  4.5  where 
we  compare  synthetics  from  several  crustal  models. 

4.3  Lg  IN  THE  TIBETAN  PLATEAU 

The  Tibetan  Plateau  is  an  especially  interesting  tectonic 
region  for  many  reasons,  but  the  one  that  concerns  us  is  that  Lg  is 
not  observed  for  paths  crossing  this  region.  This  is  pointed  out  by 
Ruzaikin,  et  al_.  (1977),  who  also  summarize  several  arguments  for 
why  the  Lg  is  absent. 

These  fall  into  two  classes: 

1.  The  waveguide  for  Lg  is  disrupted  at  the 
margins  of  the  Plateau,  or  does  not  exist  there 
at  all . 

2.  The  attenuation  in  the  crust  is  very  high 
beneath  Tibet. 

To  investigate  Lg  propagation  in  this  region,  we  compute 
theoretical  seismograms  for  two  models  for  the  Tibetan  Plateau. 
Both  are  taken  from  Chun  and  Yoshii  (1977)  and  are  listed  in  Table 
7.  The  first  model  ( TP  1 )  has  been  slightly  modified  by  adding  1  km 
of  lower  velocity  material  at  the  surface. 

The  two  Tibetan  Plateau  models  in  Table  7  both  have  a  68  km 
crust  and  the  same  half space  representing  the  upper  mantle.  The 
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Figure  16.  Synthetic  seismograms  are  shown  for  several  depths  and  ranges  in  the  oceanic  crustal 
model  of  Table  6.  The  source  is  a  strike-slip  double-couple  with  a  step  function 
time  history  with  Mo  =  10^5  dyne-cm.  The  WllSSfl  short  period  seismometer  response 
is  included.  The  source  depth  is  measured  from  the  ocean  bottom. 


TABLE  7 

TIBETAN  PLATEAU  MODELS  FROM  CHUN  AND  YOSHII  (1977) 


Layer  Depth  Thickness  a  B  p  Q 

(km)  (km)  (km/sec)  (km/sec)  (gm/cm3) 


Tibetan 

Plateau 

Model  TP1 

1 

1.0 

1.0 

3.70 

2.16 

2.10 

50 

2 

3.5 

2.5 

4.50 

2.60 

2.40 

100 

3 

13.5 

10.0 

5.98 

3.45 

2.80 

250 

4 

38.0 

24.5 

5.98 

3.45 

2.80 

1200 

5 

68.0 

30.0 

6.30 

3.64 

2.90 

2000 

6 

7.70 

4.45 

3.30 

3000 

Tibetan 

Plateau  Model  TP3 

1 

1.0 

1.0 

4.50 

2.60 

2.40 

40 

2 

3.5 

2.50 

4.50 

2.60 

2.40 

150 

3 

28.0 

24.50 

5.98 

3.42 

2.80 

200 

4 

38.0 

10.0 

5.80 

3.35 

2.75 

200 

5 

68.0 

30.0 

6.30 

3.64 

2.90 

1000 

6 

7.70 

4.45 

3.30 

2000 
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main  difference  is  that  TP3  has  a  low  velocity  layer  in  the 
mid-crust.  Chun  and  Yoshii  (1977)  prefer  models  with  this  feature. 
The  TP3  model  also  has  much  lower  Q  than  TP1. 

The  phase  and  group  velocity  dispersion  for  the  two  models  are 
plotted  in  Figure  17.  The  number  of  modes  required  to  span  a 
particular  frequency  range  is  roughly  proportional  to  the  travel 
time  from  the  surface  to  the  top  of  the  halfspace.  For  these  models 
we  computed  50  modes,  which  is  a  complete  set  for  0  to  2  Hertz. 
Recall  that  the  oceanic  crust  required  only  46  modes  for  0  to  5 
Hertz. 

Both  models  nave  a  prominent  band  of  group  velocity  minima 
between  2.9  and  3.5  Km/sec.  This  is  precisely  the  characteristic 
that  we  associate  wi;h  a  prominent  Lg  phase. 

Syntnetic  seismograms  for  the  two  models  are  shown  in  Figure 
13.  The  same  source  used  for  the  oceanic  model  synthetics  in  Sec¬ 
tion  3.2  was  used.  This  is  a  strike-slip  double-couple  with  a  step 

25 

function  time  history  with  a  moment  of  10  dyne-cm.  The  seismo¬ 
grams  for  the  two  models  are  quite  similar.  The  main  difference  is 
in  tne  amplitude,  with  the  lower  Q  TP3  having  amplitudes  that  are 
aoouc  a  factor  of  5  larger  than  those  for  TP1,  and  a  corresponding 
smft  to  lower  frequencies. 

Tne  synthetic  seismograms  have  an  Lg  arrival  that  becomes  more 
prominent  with  increasing  depth.  The  best  Lg  is  for  the  50  km 
source  depth.  At  the  shallower  depths,  the  Lg  is  rather  emergent, 
but  reaches  its  largest  amplitude  before  3.4  km/sec.  Seismograms  in 
these  models  will  be  compared  to  those  from  other  models  in  Section 
4.5 


4.4  Lg  IN  A  MODEL  FOR  THE  CENTRAL  U.S. 

4.4.1  Introduction 

In  our  semi-annual  report  (Bache  et  al . ,  1980)  we  were 
concerned  with  parametric  variations  of  a  model  intended  to 
represent  the  eastern  U.S.  The  variations  were  primarily  in  the  Q 
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Model  TP1 


Model  TP  3 


Figure  17.  Phase  and  group  velocity  dispersion  for  fifty 
modes  in  two  Tibetan  Plateau  crustal  models. 


MODEL  DEPTH  (km) 


Synthetic  seismograms  are  shown  (R  =  100  km)  for  the  two  Tibetan  Plateau  crustal 
models  of  Table  7.  The  source  is  a  strike-slip  double-couple  (Mn  =  10^  dyne-cm 
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structure  (only  frequency-independent  Q  was  considered),  but  some 
velocity  structure  variations  were  also  considered.  However,  all 
the  models  studied  had  the  same  crustal  thickness  (34  km)  and 
ha  If space  representing  the  upper  mantle.  There  are  other  velocity 
models  that  may  be  more  appropriate  for  the  eastern  and  central 
U.S.,  and  in  this  section  we  will  examine  the  Lg  predicted  for  such 
a  model.  Comparing  among  different  crustal  models  also  gives  some 
idea  of  the  dependence  of  Lg  on  the  details  of  the  velocity 
structure. 

4.4.2  Effect  of  Crustal  Thickness 

The  central  U.S.  crustal  models  to  be  considered  are  listed  in 
Table  8.  First  we  have  the  model  SI  from  our  semi-annual  report. 
It  has  a  34  km  crustal  thickness.  The  second  model  is  the  model  fit 
to  Rayleigh  waves  by  McEvilly  (1964).  This  model  has  a  38  km  crust 
and  has  been  modified  by  adding  some  lower  velocity  layers  in  the 
top  4.1  km.  The  Q  model  is  nearly  the  same  as  that  for  SI.  We  will 
call  this  model  Ml.  The  velocity  structure  is  that  used  in  our 
discussion  of  frequency-dependent  Q  models  in  Section  III. 

The  phase  and  group  velocity  dispersion  for  these  two  models 
are  plotted  in  Figure  19.  They  are  very  similar,  with  the  main 
difference  being  that  the  stationary  portions  of  the  group  velocity 
curves  are  at  somewhat  higher  velocities  for  the  McEvilly  model  than 
for  SI.  The  band  of  stationary  phases  near  3.5  km/sec  is  associated 
with  Lg.  For  SI  this  band  is  between  2.75  and  3.5  km/sec.  The 
comparable  band  for  Ml  is  between  2.9  and  3.6  km/sec.  This  leads  us 
to  expect  the  Lg  to  have  an  earlier  arrival  time  for  the  latter 
model . 

Synthetic  seismograms  are  compared  for  the  two  models  in 
Figure  20.  As  expected,  the  Ml  model  seismograms  have  an  "Lg" 
arrival  at  an  earlier  time  than  the  seismograms  for  SI.  The 
amplitudes  are  also  interesting.  At  5  and  25  km  depths  the  SI  and 
Ml  seismograms  have  nearly  the  same  maximum  amplitudes.  However,  at 
the  15  km  depth  the  maximum  Ml  amplitude  is  only  half  that  for  SI. 
We  will  plot  the  "Lg"  amplitudes  for  these  seismograms  and  two  other 
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TABLE  8 


CENTRAL  AND  EASTERN  UNITED  STATES  CRUSTAL  MODELS 


Layer 

Depth 

(km) 

Thickness  a 

(km)  (km/sec) 

(km/sec)  (gm/cm^) 

Q 

(Model  SI 

Bache 

,  Swanger  and  Shkoller, 

1980) 

1 

0.6 

0.6 

3.70 

2.16 

2.10 

20 

2 

2.6 

2.0 

4.55 

2.54 

2.20 

50 

3 

4.1 

1.5 

5.6 

3.14 

2.65 

250 

4 

6.2 

2.1 

6.1 

3.30 

2.85 

400 

5 

13.2 

7.0 

6.3 

3.41 

2.94 

1200 

6 

19.0 

5.8 

6.4 

3.46 

3.00 

1500 

7 

34.0 

15.0 

6.6 

3.59 

3.05 

2000 

8 

8.1 

4.52 

3.35 

2000 

Anisotropic  Central 

United  States  Model  (McEvilly, 

1964) 

1 

Q.6 

0.6 

3.70 

2.16 

2.10 

20 

2 

2.6 

2.0 

4.55 

2.54 

2.20 

50 

3 

4.1 

1.5 

5.60 

3.14 

2.65 

250 

4 

11.0 

6.9 

6.10 

3.50 

2.70 

4  00 

5 

20.0 

9.0 

6.40 

3.68 

2.90 

1400 

6 

38.0 

18.0 

6.70 

3.67 

2.90 

2000 

7 

00 

OC 

8.15 

4.67 

3.30 

2000 

I 


models  for  tne  crust  in  the  central  and  eastern  U.S.  at  the  end  of 
this  section. 

4.4.3  Effect  of  the  Shear  Velocity  at  the  Base  of  the  Crust 

Another  model  that  is  closely  related  to  Ml  was  discussed  in 
Section  3.3.3.  This  is  the  model  obtained  by  McEvilly  (1964)  to 
simultaneously  fit  Love  and  Rayleigh  wave  data,  rather  than  Rayleigh 
wave  data  which  were  used  to  infer  Ml.  For  the  crust,  the  only 
difference  is  that  the  3.67  km/sec  layer  above  the  mantle  is 

replaced  by  a  3.94  km/sec  layer.  Leaving  the  Q  unchanged  from  Ml, 
we  call  the  new  model  M2. 

The  phase  and  group  velocity  dispersion  for  models  Ml  and  M2 
are  compared  in  Figure  21.  As  was  mentioned  before,  the  band  of 
stationary  group  velocities  is  between  2.9  and  3.6  km/sec  for  Ml. 
Increasing  the  shear  velocity  at  the  base  of  the  crust  (Model  M2) 
causes  this  band  to  move  to  3.0  to  3.8  km/sec. 

Seismograms  for  the  model  M2  are  plotted  in  Figure  22.  As 
expected,  the  energy  is  shifted  to  much  earlier  group  arrival 

times.  Seismograms  for  the  model  M2  are  plotted  in  Figure  22.  As 
expected,  the  energy  is  shifted  to  much  earlier  group  arrivals. 

The  energy  distribution  displays  presented  in  the  previous 
report  (Bache  et  al . ,  1980)  suggested  that  nearly  all  of  the  energy 
in  the  Lg  phase  is  confined  to  the  one  or  two  deepest  layers  of  the 
crust.  This  study  suggests  that  there  is  a  direct  link  between  the 
Lg  onset  velocity  and  the  shear  velocity  at  the  base  of  the  crust. 
Examination  of  the  synthetic  Lg  onset  times  for  all  crust  models 

used  in  this  study  and  the  previous  study  suggest  an  Lg  onset 

velocity  between  96  to  98  percent  of  the  shear  velocity  in  the  lower 
crust.  The  low  value  (96  percent)  was  obtained  for  the  Tibetan 
plateau  for  which  the  synthetic  Lg  is  rather  emergent  and  onset  time 
is  not  wel 1  measured. 

It  appears  that  when  Lg  has  a  sharp  onset  the  velocity  at 
onset  may  be  very  useful  for  inferring  the  shear  velocity  of  the 
low^r  crust  over  the  path  of  propagation. 
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Model  Ml 


Model  M2 


gure  ti.  Synthetic  seismograms  are  plotted  for  the  strike-slip  doubl  e-couple  source  in 
model  M2. 


4.4.4  Lg  in  a  Central  U.S.  Model  Derived  from  Body  Wave  Refraction 

Data  ~ 

Our  attention  has  been  directed  mainly  to  models  derived  from 
observations  of  surface  waves,  usually  the  fundamental  mode.  We  now 
consider  a  model  suggested  by  Herrmann  and  Fischer  (1978)  for  the 
Central  U.S.  This  model  was  derived  to  fit  some  refraction  data 
from  an  area  near  the  Ill inoi s-Missouri  border.  Only  arrivals  with 
phase  velocities  greater  than  5.5  km/sec  were  included,  so  it  is  a  P 
wave  model.  It  is  interesting  to  examine  the  characteristics  of  Lg 
in  this  model. 

The  “Modified  Central  U.S."  model  of  Herrmann  and  Fischer 
(1978)  is  listed  in  Table  9,  after  altering  top  4.1  km  to  be  the 
same  as  SI,  Ml  and  M2.  The  Q  model  is  also  consistent  with  that  for 
these  earlier  models.  The  main  unique  feature  of  the  model  in  Table 
9,  which  we  call  Cl,  is  that  there  is  a  gradual  transition  at  layer 
boundaries. 

The  phase  and  group  velocity  dispersion  for  the  model  in  Table 
9  are  plotted  in  Figure  23.  Much  like  M2,  the  stationary  group 
velocities  are  clustered  at  high  velocities,  in  this  case  in  the 
band  between  3.1  and  3.7  km/sec.  The  synthetic  seismograms  for  this 
source  are  plotted  in  Figure  24.  The  Lg  onset  occurs  at  about  3.75 
km/sec. 

4.4.5  Amplitudes  of  the  Synthetic  Lg  for  the  Central  and  Eastern  U.S. 

For  each  of  the  four  models  discussed  in  this  section,  we  have 
computed  synthetic  seismograms  with  the  same  source.  The  Q  models 
are  also  very  much  the  same  for  each  model;  the  main  differences  are 
in  the  velocity  structure.  On  each  seismogram  we  measured  the  Lg 
amplitude  which  is  taken  to  be  the  maximum  sustained  amplitude 
within  10  seconds  of  the  3.5  km/sec  group  arrival  time.  These 
amplitudes  are  plotted  in  Figure  25. 

4.5  COMPARISON  OF  SYNTHETIC  Lg  SEISMOGRAMS  IN  DIFFERENT  CRUSTAL 

STRUCTURES 

We  have  been  computing  synthetic  Lg  seismograms  in  three 
distinctly  different  tectonic  regions.  In  Section  4.2,  it  was  the 


TABLE  9 


Cl  8ASE  ON  THE  MODIFIED  CENTRAL  U.S.  MODEL 
(HERRMANN  AND  FISCHER,  1978) 


Depth  Thickness  a  s  p 


Layer 

(km) 

(km) 

(km/sec) 

(km/sec)  (gm/cm3) 

Q 

1 

0.6 

0.6 

3.70 

2.16 

2.10 

20 

2 

2.6 

2.0 

4.55 

2.54 

2.20 

50 

3 

4.1 

1.5 

5.60 

3.14 

2.65 

250 

4 

7.0 

2.9 

6.10 

3.52 

2.56 

400 

5 

8.0 

1.0 

6.18 

3.57 

2.59 

400 

6 

9.0 

1.0 

6.32 

3.65 

2.65 

400 

7 

18.0 

9.0 

6.40 

3.70 

2.68 

1400 

8 

19.0 

1.0 

6.48 

3.74 

2.71 

1600 

9 

20.0 

1.0 

6.62 

3.82 

2.76 

1800 

10 

38.0 

18.0 

6.70 

3.87 

2.79 

2000 

11 

39.0 

1.0 

7.06 

4.08 

2.93 

1500 

12 

40.0 

1.0 

7.79 

4.50 

3.20 

1800 

13 

oo 

00 

8.15 

4.71 

3.34 

2000 
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Figure  25. 
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deep  ocean  and  in  Section  4.3  it  was  the  Tibetan  Plateau  where  the 
crust  is  very  thick.  In  Section  4.4  we  studied  several  models  for 
the  eastern  U.S.  wnich  represents  a  stable  platform  area.  In  this 
section  we  directly  compare  seismograms  for  these  regions. 

The  first  comparisons  are  shown  in  Figure  26.  The  seismograms 
were  computed  for  the  standard  strike-slip  double-couple  source  we 
have  been  using,  and  at  two  source  depths.  The  oceanic  crustal 
model  is  very  much  different  from  the  others  in  having  the  energy 
distributed  in  a  wide  group  velocity  window.  These  is  clearly  no 
waveguide  for  Lg  in  the  oceanic  crust,  and  this  is  why  it  does  not 
propagate  there. 

The  Tibetan  Plateau  model  gives  interesting  results.  For  the 
shallow  source  it  has  an  Lg  phase  just  as  prominent,  if  not  more  so, 
that  the  stable  platform  models.  However,  with  a  source  depth  of  5 
km,  the  Lg  phase  nearly  disappears  in  the  Tibetan  Plateau  model.  We 
need  to  compare  the  seismograms  at  greater  ranges  to  draw  firm 
conclusions  about  the  TP3  model  as  a  propagator  of  Lg. 

In  Figure  27,  we  compare  synthetic  seismograms  at  three  source 

depths  in  five  crustal  models.  These  are  the  Tibetan  Plateau  Model 

TP3  and  four  different  models  representing  the  eastern  and  central 

U.S.  The  source  is  a  strike-slip  double-couple  with  a  step-time 
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history  with  moment  10  dyne-cm.  The  models  are  described  in 
Sections  4.3  and  4.4. 

In  summary,  Lg  fails  to  propagate  in  oceanic  structure  because 
there  is  no  waveguide  for  the  phase  quite  like  that  in  continental 
crusts.  For  the  particular  crustal  models  used  for  our  computa¬ 
tions,  energy  of  the  multi-mode  Rayleigh  waves  a^e  distributed  over 
a  much  wider  group  velocity  window  than  is  normally  associated  with 
Lg.  The  large  dispersion  also  results  in  considerably  lower 
amplitudes  in  the  group  velocity  window  of  interest  than  those  in 
continental  crusts.  The  major  structural  feature  lacking  in  the 
oceanic  model  used  is  a  distinct  Moho.  There  are  particular 
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synthetic  seismograms  (R  =  200  km)  for  four  crustal  models. 
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Figure  27.  Synthetic  seismograms  are  compared  at  three  source  depths  in  a  Tibetan  Plateau 
crustal  model  (TP3)  and  four  models  for  the  central  and  eastern  United  States. 
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locations  where  a  distinct  Moho  appears  to  be  present  (D.  V. 
Helmoerger,  personal  communication).  At  such  places  it  may  be 
possible  to  generate  an  Lg-like  phase. 


Standard  Tibetan  Plateau  models  provide  an  excellent  waveguide 
for  Lg.  The  reason  it  does  not  propagate  in  this  region  must  be 
sought  elsewhere.  One  other  possible  explanation  is  that  the  Lg 
propagates  in  the  region,  but  is  severely  attenuated  upon  passage 
across  the  boundary  with  other  regions  where  attempts  to  observe  Lg 
are  made.  Another  explanation  is  that  there  is  an  unusual  Q 
structure  that  attenuates  this  phase. 

For  all  crustal  models  considered  here  in  which  synthetic  Lgs 
are  prominent,  the  onset  group  velocity  of  the  phase  is  about  96-98 
percent  of  the  shear  velocity  at  the  base  of  the  onset.  The  close 
ties  between  Lg  onset  velocities  and  lower  crust  shear  velocities  is 
to  be  expected  given  the  energy  distributions  of  Lg  found  in  Bache 
et  al ,  (1980).  When  one  considers  what  Lg  primarily  is,  namely, 
critically  reflected  shear  waves  guided  between  the  Moho  and  earth's 
surface,  one  can  make  a  more  precise  statement  about  the  onset 
time.  The  onset  time  will  be  that  of  the  arrival  of  the  primarily 
snear  wave  wide-angle  reflection  off  the  Moho.  Because  of  crustal 
velocity  distribution  with  depth,  this  time  is  very  nearly  that  of  a 
snear  wave  propagating  horizontally  through  the  lower  crustal 
layering. 


V.  DEPENDENCE  OF  Lg  ON  SOURCE  DEPTH 


In  our  semi-annual  report  (8ache  et  al . t  1980)  we  discussed 
the  theoretical  dependence  of  Lg  on  source  depth,  using  synthetic 
seismograms  in  the  eastern  U.S.,  crustal  model  SI  (Table  8).  Those 
results  were  not  entirely  satisfactory,  because  a  complete  and 
consistent  set  of  synthetic  seismograms  was  not  available  to 
construct  amplitude-depth  curves.  Since  that  time,  we  have  computed 
a  consistent  set  and  the  new  results  will  be  presented  here. 

Synthetic  Lg  seismograms  were  computed  for  a  series  of  source 
depths  in  the  model  SI.  The  source  was  a  point  double-couple  at  the 
three  fundamental  orientation:  strike-slip,  vertical  dip-slip  and 

45°  dip-slip.  The  source  time  history  was  taken  to  be  a  step 

25 

function  with  a  fixed  moment  of  10  dyne-cm.  The  WWSSN  short 
period  seismometer  response  was  included.  Synthetics  were  computed 
at  a  range  of  1000  km  for  all  three  source  orientations.  In 
addition,  the  calculations  were  done  at  500  km  for  the  strike-slip 
douole-couple. 

The  synthetic  seismograms  are  shown  in  Figure  28.  The  Lg 
amplitude  was  measured  on  each  seismogram.  This  was  taken  to  be  the 
maximum  sustained  amplitude  (at  least  3  cycles)  within  ten  seconds 
of  the  3.5  km/sec  group  arrival  time. 

The  amplitudes  are  plotted  versus  source  depth  in  Figure  29. 
This  is  basically  a  cleaner  version  of  a  similar  plot  that  was  shown 
in  our  semi-annual  report.  The  computations  are  different  in  that: 

(1)  More  depths  are  sampled, 

(2)  All  synthetics  are  computed  with  the  same 
instrument  response, 

(3)  A  more  consistent  and  conventional  technique 
for  measuring  Lg  amplitude  has  been  employed. 

While  the  results  are  different  in  detail,  the  trends  are  the 
same  as  shown  in  our  semi-annual  report.  The  strike-slip  and  45° 
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Figure  28.  (continued) 
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Figure  28.  (continued) 
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Figure  28.  (continued) 


Vertical  Strike-Sli 


Lg  amplitudes  as  a  function  of  source  depth  for  fixed  seismic 
moment  at  1000  km  in  an  EUS  crustal  model. 


dip-slip  orientations  show  a  similar  decrease  of  Lg  amplitude  with 
depth.  Most  of  the  decrease  is  at  the  shallow  and  deep  ends  of 
depth  interval.  Between  3  and  21.5  km  depths  the  decrease  is  only  a 
factor  of  two  or  so.  The  normal  dip-slip  fault  exhibits  the 
opposite  trend,  except  near  the  surface. 

With  this  depth  sampling,  we  can  see  that  the  proximity  to  an 
interface  might  be  making  some  difference.  The  somewhat  anomalous 
values  at  depths  of  9.5  and  18  km  are  probably  influenced  by  their 
distance  from  the  interfaces. 

The  range  dependence  of  the  depth  effect  is  indicated  by 
comparing  the  R  =  1000  km  and  R  =  500  km  cases  for  the  strike-slip 
orientation.  To  parameterize  this  effect,  we  compute  the  effective 
attenuation  (rSyn)  'for  each  depth,  using  the  procedure  introduced 
by  Nuttli  (1973)  (see  Section  3.2).  These  are: 


Depth 

Ysynfdeg"1) 

1 

0.19 

3 

0.28 

5 

0.28 

7 

0.34 

9.5 

0.30 

12 

0.34 

15 

0.37 

18 

0.21 

21.5 

0.26 

25 

0.41 

There  are  no  striking  trends  in  these  y*  They  have  a  mean  of  0.30 
deg”*  (0.07  standard  deviation).  This  is  close  to  the  value  we 
found  for  the  structure  SI  with  an  explosion  source. 

The  amplitudes  in  Figure  29  are  for  a  constant  Mq  source. 
The  relationship  between  Mq  and  is  independent  of  the  elastic 
properties  of  the  source  layer,  so  the  plot  is  also  for  a  constant 
Ms  source.  It  may  also  be  of  interest  to  plot  the  depth 
dependence  for  a  constant  m^.  To  show  this,  it  is  necessary  to 
scale  by  the  elastic  properties  of  the  source  layer. 


Using  geometrical  ray  theory  for  a  spherical  earth  with  depth 
dependent  velocity,  the  far-field  P-wave  amplitudes,  ignoring  source 
radiation  pattern,  can  be  written  ( Aki  and  Richards,  1980,  Chapter  9) 


Up(rR,  a,  t)  - 


cos  iR  MQ(t  -  T) 


4ir 


- 5- 

PR  ps  °R  °s 
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R(  Xr, 


Is) 


wi  th 


R(rR,  r$)  = 


rR  rs 


cos  iR  cos  is  sin  al-j—j/P 


1/2 


where 

M  =  moment 
0 

T  *  travel  time 

i  =  takeoff  or  incident  angle  from  the  vertical 

r  =  radius  from  earth's  center 

a  =  range  (angle) 

P  =  ray  parameter 

and  S,  R  refer  to  source  and  receiver. 

The  factor  depending  on  source  properties  is 


1 _ 

cos  i 


s 
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For  crustal  source  depths  and  teleseismic  distances,  r$  and 
cos  i  vary  by  only  a  few  percent  so  the  dominant  depth  effect  is 
the  factor 
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p  a. 


The  computed  Ig  amplitudes  for  fixed  mb  are  shown  in 
Figure  30.  Note  that  much  of  the  depth  dependence  for  vertical 
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Oblique  Thrust  t  ^ Vertical  Str 


Figure  30.  Lg  amplitudes  as  a  function  of  source  depth  for  fixed  "iti| 
1000  km  ir  an  EDS  crustal  model. 


strike-slip  and  oblique  thrust  sources  is  compensated  for  by  the 
depth  dependent  factor  above.  For  vertical  dip-slip  sources,  on  the 
other  hand,  the  depth  effect  becomes  even  more  pronounced. 
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VI.  THREE-DIMENSIONAL  SOURCE-SIMULATION  CAPABILITIES 


Some  important  seismic  source  phenomena  are  fundamentally  both 

three-dimensional  and  nonlinear.  Tectonic  earthquakes  are  an 

example.  Buried  explosions  near  nonhorizontal  geologic  interfaces 

or  near  imperfectly-welded  rock  joints  represent  another  class  of 

three-dimensional  sources.  One  tool  for  studying  the  near-  and 

far-field  ground  motion  excited  by  such  sources  is  three-dimensional 

finite  difference  simulation.  For  this  purpose,  Systems,  Science 

and  Software  (S3)  developed  the  TRES  code  (Cherry,  1977),  a 

3 

small -strain  finite  difference  code  which  operates  on  the  S 
UNIVAC  1100/81.  TRES  was  constructed  to  optimally  exploit  the 
scalar  processing  capabilities  of  the  UNIVAC. 

Because  of  the  large  computational  requirements  of  three-dimen¬ 
sional  seismic  modeling,  the  TRES  algorithm  was  subsequently 

reprogrammed  (by  the  Institute  for  Advanced  Computations)  to  run  on 
the  ILLIAC  IV  parallel -processing  computer.  This  conversion  was  a 
multi -man-year  effort,  because  of  the  unique  demands  of  the  ILLIAC. 
Code  development  for  the  ILLIAC  proved  to  be  about  an  order  of 
magnitude  more  time-consuming  than  comparable  development  efforts  on 
the  UNIVAC.  Once  in  place,  however,  the  ILLIAC  version  of  TRES, 
designated  I4TRES,  executed  three-dimensional  dynamic  simulations, 
involving  up  to  about  one  million  finite  difference  nodes,  in  a  few 
hours  of  computing  time.  This  code  was  used  by  S  during  FY  '79 

to  develop  a  shear  failure  model  for  earthquakes,  and  during  FY  '80 

to  perform  several  large-scale  earthquake  simulations,  with  the 
objective  of  improving  the  theoretical  basis  for  seismic  source 
identification  techniques. 

Two  tasks  under  this  contract,  allocated  approximately  ten 
man-weeks  each,  were  directed  toward  improving  this 
three-dimensional  source-modeling  capability.  Task  1  was  to  provide 
support  to  Technology  Development  Corporation  (TDC),  (formerly  the 
Institute  for  Advanced  Computation)  in  order  to  modify  the  I4TRES 
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code.  Task  2  involved  converting  TRES  from  the  S  UNIVAC  to  the 
CRAY  I  computer  at  the  National  Center  for  Atmospheric  Research 
(NCAR),  and  then  comparing  the  performance  of  TRES  on  the  UNIVAC, 
CRAY,  and  ILLIAC  computers.  Below,  we  summarize  the  work 
accomplished  under  these  tasks  ind  our  conclusions  about  the 
three-dimensional  modeling  capabilities  of  the  CRAY  I  version  of 
TRES. 

6.1  SUMMARY  OF  TASK  1 

3 

Under  Task  1,  S  provided  TDC  with  appropriate  algorithms 
and  test  problems  designed  to  enable  TDC  to  upgrade  the  modeling 
capability  of  I4TRES.  The  objectives  were: 

1.  Complete  and  test  the  multiple-materials  version  of  TRES 
which  had  been  under  development  by  TDC. 

2.  Make  it  compatible  with  S^'s  nonlinear  rupture  model. 

3.  Allow  for  general  free-surface  interactions,  including 
outcropping  thrust  faulting. 

3 

The  S  algorithms  and  test  problems  are  documented  in  Appendixes  A 
and  B. 


As  of  the  begining  of  January,  1981,  Objective  1  had  been 
successfully  met  as  evidenced  by  the  results  of  two  key  test 

problems  discussed  below.  Objective  2  was  successfully  met  in 
March,  1981.  Completion  of  Objective  3  has  not  been  successfully 

demonstrated. 

I4TRES  was  drastically  reformulated  by  TDC  to  accommodate 
material  heterogeneity.  Test  Problem  I,  described  in  Appendix  B, 
was  intended  to  verify  that  the  revised  code  correctly  modeled 

faulting  in  a  homogeneous,  elastic  half space.  In  particular,  we 
wanted  to  verify  that  the  revised  difference  equations  reduce  to  the 
correct  boundary  conditions  at  the  free  surface  and  at  the  fault 

plane. 
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Figure  31  shows  free-surface  displacement  time-histories 
computed  by  I4TRES  for  Problem  I.  Also  shown  are  two  independent 
solutions,  obtained  by  the  Cagniard-de  Hoop  method  (Johnson,  1974) 
and  a  quasi-three-dimensional  finite  element  method  (Day,  1977), 
respectively.  Comparison  of  the  Cagniard-de  Hoop  and  finite  element 
solutions  illustrates  the  filtering  effect  associated  with  discrete 
methods  such  as  finite  element  and  finite  difference:  frequencies 
corresponding  to  shear-wavelengths  shorter  than  about  six 
zone-dimensions  are  inaccurately  propagated  by  the  numerical  methods 
and  have  been  suppressed  by  artifical  viscosity.  Since  the  zone 
size  used  in  I4TRES  was  twice  as  large  as  that  used  in  the  finite 
element  solution,  the  filtering  effect  is  even  more  severe  for  the 
I4TRES  solution,  as  expected.  However,  low-frequency  character¬ 
istics  are  preserved  in  the  I4TRES  solution,  and  demonstrate  that 
the  code  is  working  satisfactorily  for  this  test  case. 

Test  Problem  II  was  designed  to  verify  that  the  revised  I4TRES 
correctly  models  wave  propagation  in  a  non-homogeneous  halfspace. 
Even  for  the  case  where  the  halfspace  consists  of  a  few  discrete, 
horizontal,  elastic  layers,  only  a  few  complete  solutions,  valid  in 
the  near  field,  are  available  in  the  literature.  Problem  II  treats 
a  problem  for  which  we  have  a  published  solution  obtained  by  the 
quasi-three-dimensional  finite  element  method  (Day,  1977).  Again, 
the  three-dimensional  I4TRES  solution  and  quasi -three-dimensional 
finite  element  solution  must  be  compared  in  light  of  the  fact  that 
the  zone  size  of  the  I4TRES  calculation  was  twice  that  of  the  finite 
element  calculation.  Thus,  the  I4TRES  solution  should  be 
essentially  a  low-pass  filtered  version  of  the  finite  element 
solution.  Figure  32  shows  that  this  is  the  case. 

Test  Problem  III  was  designed  to  ensure  that  the  nonlinear 
rupture  model  operates  correctly  in  conjunction  with  the 
multiple-material  version  of  I4TRES.  Successful  completion  of  this 
test  problem  was  verified  by  comparing  the  new  numerical  results  to 
numerical  results  obtained  for  the  same  problem  with  an  earlier 
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Figure  31.  Free  surface  displacement  waveforms  for  Test  Problem  I,  computed 

with  I4TRES  (the  ILLIAC  IV  version  of  the  TRES  code);  displacements 
computed  by  two  alternate  methods  are  shown  for  comparison.  The 
slip  functions  shown  were  applied  to  a  1  km?  fault  area,  centered 
at  5  km  depth.  Observers  are  at  an  azimuth  of  22. 5"-  from  epicentral 
distance. 
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Figure  32.  Free  surface  displacements  for  I4TRES  Test  Problem  II,  compared 
to  a  two-dimensional  finite  element  solution  (Day,  1977).  Source 
receiver  conf iguration  is  the  same  as  for  Figure  31. 


version  of  I4TRES.  Discrepancies  were  found  to  be  insignificant  and 
fully  attributable  to  the  (algebraically  equivalent)  rearrangement 
of  the  difference  equations  in  the  new  version  of  I4TRES. 

In  conclusion,  the  results  of  Test  Problems  I,  II  and  III  are 

very  satisfactory.  The  I4TRES  code  as  now  configured  can  accurately 

treat  three-dimensional  wave-propagation  in  a  heterogeneous  medium, 

including  the  stress-relaxation  which  accompanies  fault  propagation 

3 

at  a  specified  rupture  velocity.  In  addition,  the  S  shear 
failure  algorithm  given  in  Appendix  A  has  been  successfully 
integrated  into  the  I4TRES  code.  The  capability  to  model 
non-vertical  faults  has  not  been  demonstrated  by  successful 
completion  of  Test  Problem  IV. 

6.2  SUMMARY  OF  TASK  2 

The  objective  of  Task  2  was  to  evaluate  the  performance  of  the 
TRES  code  on  the  UNIVAC  1100/81,  ILLIAC  IV,  and  CRAY  I  computers. 
The  main  accomplishments  were: 

1.  Creation  of  a  version  of  the  TRES  code  in  CRAY  FORTRAN 
(CFT) . 

2.  Partial  vectorization  of  TRES  to  exploit  the 
vector-processing  capabilities  of  the  CRAY  I. 

3.  Comparison  of  TRES  code  execution  times  on  the  three 
computers. 

Conversion  of  the  vectorization  of  TRES  are  briefly  discussed 
in  Appendix  C,  The  remainder  of  this  section  gives  the  results  of 
the  comparison  of  TRES  runs  on  the  three  computers. 

It  is  difficult  to  make  precise  comparisons  of  run  times  among 
the  three  computers.  The  real  objective  of  such  a  comparison  should 
be  to  try  to  predict  the  relative  performances  of  the  computers  for 
large  "production"  runs.  For  this  purpose,  the  most  useful  form  in 
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which  to  express  the  timings  would  be  as,  say,  CPU-seconds  per 
finite-difference  zone  per  time  step.  In  fact,  once  this  figure  is 
deduced  for  a  mode rate- si  zed  test  problem,  we  can  predict  UNIVAC 
run-times  with  reasonable  confidence  for  very  large  problems  simply 
by  assuming  that  run-time  is  proportional  to  the  product  of  the 
number  of  zones  and  the  number  of  time  steps  in  the  problem. 

However,  computing  times  for  the  ILL IAC  and  CRAY  machines  are 
distinctly  nonlinear  in  the  number  of  zone-cycles  and  the  form  of 
the  nonlinearity  is  different  for  each.  Furthermore,  for  both  the 
UNIVAC  and  CRAY  computers,  it  is  essential  to  correct  test  problem 
timings  for  the  fact  that  only  part  of  the  finite  difference  grid  is 
active  during  the  initial  stages  of  a  given  calculation. 

Because  of  the  above  considerations,  we  determined  that  a 
single  test  problem  was  inadequate  to  enable  us  to  predict  the 
relative  performance  of  the  computers.  Instead,  we  performed 
several  calculations  with  various  grid  configurations.  Based  on 
these  runs,  we  compare  the  expected  performance  of  the  three 
computers  for  problem  sizes  of  practical  interest.  Table  10 
summarizes  the  runs  used  in  the  analysis. 

6.2.1  ILLIAC  Computing  time 

TRES  run  times  on  this  machine  are  reasonably  well  understood 
as  a  result  of  several  years  experience  with  large-scale 
calculations.  Let  J,  K,  and  L  be  the  number  of  finite  difference 
zones  in  the  x,y,  and  z  directions  so  that  the  total  number  of  zones 
in  the  problem  is  the  product  JKL.  Then  the  computer  time  (in 
seconds)  per  cycle,  T,  is  approximately: 


T 


.094  L 


K  1 

J 

[m 

(I) 


where  the  bracket  symbol  [x]  indicates  the  smallest  integer  greater 
than  or  equal  to  x.  Thus,  T  is  proportional  to  L  and  independent  of 
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Is  the  average  fraction  of  the  grid  estimated  to  be  active  during  a  typical  production  run. 
C  =  I  0.75  for  UN  I VAC  and  CRAY 


J  up  to  J  =  60.  For  very  small  problems  (k  <  20),  T  is  also 
independent  of  K;  while  for  very  large  problems  (k  »  20),  T  is 
practically  proportional  to  K.  This  J  and  K  dependence  renders  the 
ILLIAC  relatively  much  more  efficient  for  very  large  problems. 
Combined  hardware/software  restrictions  limit  ILLIAC  problem  size  to 
J  =  80,  K  =  80,  L  =  160. 

6.2.2  CRAY  Computing  Time 

Computing  time  for  the  CRAY  depends  linearly  on  K  and  L.  For 
J,  however,  dependence  is  nonlinear  due  to  the  vector  processing 
capabilites  of  the  machine  and  the  partially  vectorized 
configuration  of  the  code.  CPU-time  divided  by  J  decreases  somewhat 
with  J  up  to  J  =  64;  for  larger  J,  there  is  first  some  increase  in 
T/J,  then  a  decrease  up  to  J  =  128. 

An  activity  test  has  been  incorporated  into  TRES  so  that,  for 
most  calculations,  T  increases  in  approximate  proportion  to  the 
time-step  number  during  the  first  half  of  the  problem  (while  the 
wavefield  spreads  and  the  entire  grid  gradually  becomes  active); 
then  T  remains  stationary  at  a  maximum  value  during  the  second  half 
of  the  problem  (during  which  the  disturbance  at  the  external  grid 
boundaries  propagates  back  into  the  central  region  of  the  grid). 
Thus,  for  a  typical  problem,  3/4  of  the  grid  will  be  active  during 
tne  average  time-step.  The  effects  of  the  activity  test  must  be 
accounted  for  both  in  interpreting  the  test-problem  performance  and 
in  extrapolating  to  large-scale  problems.  After  compensating  for 
this  grid-activity  effect,  we  find  approximately 

Tasc(J)KLJ  (2) 


where 


c(60)  %  .00012 
c(21 )  =»  .00019 
c(80)  .00013 
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6.2.3 


UNIVAC  Computing  Time 

For  this  machine,  there  is  no  preferred  index;  that  is,  we 
expect  approximate  proportionality  of  T  to  JKL.  The  grid  activity 
effect  is  the  same  as  described  for  the  CRAY.  After  correcting  for 
activity,  we  find  approximately 

T%cJKL  (3) 


where 


C  as  .0019 
6.2.4  Summary 

To  compare  ILLIAC  times  to  those  of  the  other  two  machines,  we 
should  account  for  the  benefits  of  limited  grid  activity  during 
early  timesteps.  In  Figure  33,  we  plot  JKL/T,  the  number  of  zones 
times  the  number  of  timesteps  per  CPU-second.  Values  of  JKL/T 
implied  by  Equations  2  and  3  for  the  UNIVAC  and  CRAY,  respectively, 
have  been  multiplied  by  4/3  since,  as  we  have  said,  the  average 
number  of  active  zones  during  a  large  calculation  is  usually  about 
3/4  JKL. 

The  figure  indicates  that  the  CRAY  is  somewhat  faster  than  the 
ILLIAC  for  most  TRES  problems,  sometimes  by  as  much  as  a  factor  of 
two.  For  certain  optimal  conf igurations ,  however,  the  ILLIAC  can 
probably  beat  the  CRAY  by  as  much  as  15  percent  or  so. 

In  special  cases  where  the  average  active  fraction  of  the  grid 
is  nearly  1,  the  ILLIAC  and  CRAY  will  give  quite  similar  times  for 
most  configurations;  on  the  other  hand,  for  special  cases  in  which 
the  fractional  activity  is  less  than  3/4,  the  advantage  of  the  CRAY 
will  be  proportional ly  greater  than  suggested  by  Figure  33. 
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timesteps/CPU-second 


ILL IAC  (k  =  20,  40,  60,  or  80) 


I  ILIAC  (k  =  30) 


Figure  33.  Speed  of  execution  of  TRES,  in  zone-cycles  per  CPU-second.  The 
horizontal  axis,  J,  is  the  problem  size  (number  of  finite  dif¬ 
ference  zones)  in  the  X  direction.  Data  points  are  given  by 
solid  circle,  squares,  and  triangles.  Curves  represent  inter¬ 
polation  and  extrapolation  of  the  data  as  described  in  the  text 


For  large  problems,  say  J  greater  than  40,  the  CRAY  currently 
runs  IRES  approximately  15  times  as  fast  as  does  the  UNIVAC. 

At  this  point,  we  should  emphasize  once  again  that  the  TRES 
code  has  been  only  partially  vectorized  to  exploit  the  CRAY'S 
capabilites.  We  estimate  that  a  substantial  effort  to  fully 
vectorize  TRES  would  yield  another  factor  of  four  or  greater 
improvement  in  speed  on  the  CRAY. 
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VII.  SIMULATION  OF  THE  VERTICAL  GROUND  ACCELERATIONS 
IN  THE  1971  SAN  FERNANOO,  CALIFORNIA,  EARTHQUAKE 


7.1  INTRODUCTION 

In  this  section  we  summarize  an  effort  to  simulate  the 
near-field  ground  response  in  the  1971  San  Fernando,  California, 
earthquake.  A  computer  program  for  the  simulation  of  the  ground 
acceleration  provided  by  has  been  used  by  the  Mission  Research 
Corporation  (MRC)  to  estimate  the  level  of  acoustic  disturbances 
expected  for  that  event.  Simulated  ionospheric  disturbances  can  be 
compared  to  those  of  underground  tests  to  determine  the  potential  of 
ionospheric  monitoring  techniques  as  a  possible  discriminant. 

S  s  part  in  this  overall  effort  was  to  provide  an  efficient 
ground  motion  simulator  to  MRC  for  obtaining  reasonable  estimates  of 
the  near-field  motion.  Here,  we  describe  the  methods  used  in 
computer  algorithms  provided  to  W.  Wortman  of  MRC,  S.  Warshaw  of 
LLL,  and  J.  Bannister  of  Sandia  Laboratories.  In  Section  7.2,  the 
mechanism  of  the  San  Fernando  earthquake  is  briefly  reviewed,  and 
the  sections  following  thereafter  discuss  the  implementation  of  our 
algorithm.  Our  final  model  for  the  source  includes  three  distinct 
events  which  are  represented  by  three  circular  cracks  situated  on 
the  rupture  surface.  The  choice  of  the  source  parameters  of  the 
discrete  sources  is  guided  by  short  period  teleseismic  observations 
and  the  observed  near-field  velocities  at  Pacoima  Dam. 

7.2  PAST  STUDIES  OF  THE  SAN  FERNANDO  FAULTING  MECHANISM 

The  1971  San  Fernando  earthquake  is  probably  the  most  studied 
earthquake  in  history.  The  location  and  moderate  size  resulted  in 
numerous  teleseismic  and  strong-motion  data  of  reasonable  quality. 
Over  the  past  decade  numerous  studies  of  the  teleseismic  and 
near-field  data  have  been  made  using  both  short  and  long  period 
information.  For  the  present  application,  the  intermediate 
frequencies  are  of  interest.  The  acoustic  pressure  signatures  are 
most  closely  reflected  by  the  characteristics  of  near-field  ground 
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velocity  (W.  Wortman,  personal  communication) ,  which,  at  Pacoima 
Dam,  falls  into  the  frequency  range  of  about  0.5  to  3.0  Hz.  For 
this  frequency  range,  studies  of  near-field  velocity  and  short 
period  teleseismic  data  are  most  informative.  Most  notable  studies 
of  this  variety  are  Hanks  (1974),  Boore  and  Zoback  (1974),  Bache  and 
Barker  (1978),  Bouchon  (1978),  and  Cherry,  et  al .  (1976). 

From  all  of  these  studies,  there  are  a  few  fundamental 
observations  about  the  faulting  mechanism  which  can  be  identified: 

1.  The  faulting  process  was  quite  complicated.  This  was 

first  noted  by  Hanks  (1974),  who  noted  the  presence  of 
several  strong  arrivals  on  both  the  short  period 

teleseismic  body-wave  data  and  the  near-field 
accelerations. 

2.  There  was  an  unusually  high  stress  release  in  the 
vicinity  of  the  hypocenter.  Estimates  generally  range 
from  200  to  500  bars.  Consistent  findings  in  this  regard 
were  found  by  Hanks  (1974),  Bache  and  Barker  (1978)  and 
Bouchon  (1978). 

3.  The  rupture  surface  was  kinked.  This  was  apparent  when 

the  fault-plane  solution  for  the  initiation  was  found  not 

to  project  from  the  hypocenter  to  the  location  of  surface 
breakage. 

4.  There  appears  to  be  large  static  slip  at  both  the 

hypocentral  region  and  at  shallow  depths  near  the  surface 
breakage  (Alewine,  1974;  Heaton  and  Helmberger,  1978). 

5.  There  is  some  source  of  significant  high  frequency 
radiation  initiated  at,  or  near,  the  surface  breakage. 

Any  reasonable  source  model  for  the  event  must  include  most,  or  all, 
of  the  above  mentioned  features.  Features  1-3,  and  5  are  very 
important  to  determining  the  characteristics  of  the  near-field 
radiation  in  the  frequency  range  of  interest. 
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Our  interpretation  of  the  results  of  previous  studies  suggests 
a  multiple  source  mechanism  consisting  of  three  discrete  events. 
The  ground  velocities  infered  from  the  recorded  accelerograms  at 
Pacoima  Dam  indicate  three  distinct  disturbances  in  the  intermediate 
frequency  band  (Figure  34).  The  largest  velocities  are  primarily 
associated  with  the  first  phase  beginning  about  2.5  seconds  into  the 
record  on  all  three  components.  On  the  component  S74W,  two  other 
strong  phases  are  quite  obvious  at  6  and  8  seconds.  These  secondary 
arrivals  are  present  on  the  other  two  components,  but  are  less 
obvious. 

The  source  of  the  first  two  disturbances  can  be  found  in  the 
model  proposed  by  Bache  and  Barker  (1978).  Their  model,  which  was 
constrained  primarily  from  teleseismic  data  and  partially  by  the 
Pacoima  Dam  vertical  records,  consists  of  two  primary  areas  of 
stress  release  (Figure  35).  The  dominant  radiation  initiates  from 
the  hypocentral  region  with  a  large  stress  drop.  A  second  suggested 
region  of  concentrated  energy  release  is  located  just  above  the  kink 
in  the  fault.  The  arrival  times  at  Pacoima  Dam  from  this  part  of 
the  Bache  and  Barker  model  is  quite  consistent  with  the  second  burst 
of  energy  observed  6  seconds  into  the  record. 

The  third  source  of  energy  has  largely  gone  unexplained.  This 
disturbance  produced  the  large  accelerations  (over  1  g)  observed  on 
both  horizontal  components,  which  at  the  time  were  the  largest 
ground  accelerations  ever  recorded.  The  traditional  hypothesis 
concerning  the  source  of  this  radiation  is  that  it  results  from 
Rayleigh  waves  generated  by  a  surface  breakout  of  the  rupture.  The 
physics  of  such  a  process  is  not  well  understood,  and  quantitative 
verification  of  this  hypothesis  has  not  been  made.  We  feel  that 
some  treatment  of  this  disturbance  is  necessary,  but  do  not  feel 
that  the  traditional  explanation  is  necessarily  a  good  one.  First, 
the  amplitude  of  the  vertical  motion  is  quite  small  relative  to  the 
horizontal.  The  ground  velocity  on  S74W  is  roughly  three  times  that  of 
the  vertical  and  S16E  roughly  twice  the  vertical.  This  does  not 
appear  consistent  with  Rayleigh  wave  motion  in  what  is  essentially  a 
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Processed  ground  velocity  observed  at  Pacoima  Dam  (from  EERL,  1973) 


homogeneous  half-space.  Second,  static  observations  (Alewine,  1974) 
and  modeling  of  the  near-field  displacements  (Heaton  and  Helmberger, 
1978)  indicate  large  slip  (^4m)  at  depths  of  a  few  km  which  is  twice 
that  at  the  surface.  Given  these  two  observations,  we  suggest  that 
a  more  likely  source  of  the  radiation  was  a  shallow  stress 
concentration  nearly  beneath  Pacolma  Dam.  Shear  waves  radiating 
from  below  would  produce  relative  horizontal  to  vertical  motion  more 
consistent  with  the  observations.  The  center  of  this  disturbance 
could  also  correspond  to  the  location  at  maximum  slip  at  shallow 
depths. 

7.3  REPRESENTATION  OF  DISCRETE  SOURCES  OF  RADIATION 

In  this  section,  we  describe  the  source  representation 
employed  for  the  discrete  sources  of  radiation.  We  will  use  a 
modification  of  a  method  used  by  Day  et  al_;_,  (1980)  which  used  an 
approximation  to  the  far-field  radiation  from  an  expanding  circular 
shear  crack  which  stops  gradually. 

Several  closed-form  analytic  approximations  to  the  slip 
history  of  a  circular  shear  crack  with  uniform  stress-drop  are 
available.  These  are  generally  motivated  by  the  early  analytical 
solution  of  Kostrov  (1974)  and  numerical  simulations,  like  Madariaga 
(1976),  which  include  the  stopping  of  rupture.  Boatwright  (1980) 
provides  a  summary  of  available  analytic  expressions.  Even  though 
simple  expressions  for  the  slip  history  on  the  fault  exist,  simple 
representations  of  the  geometrical  far-field  radiation,  in  general, 
do  not.  The  exception  is  the  model  of  Sato  and  Hirasawa  (1973). 

The  Sato  and  Hirasawa  model  (which  we  will  call  the  S  &  H 
model)  is  simply  an  expanding,  constant  rupture  velocity,  circular 
crack  for  which  slip  everywhere  on  the  crack  terminates 
instantaneously  when  the  rupture  reaches  a  prescribed  radius.  Their 
model  is  not  a  rigorous  dynamic  solution,  as  shown  by  the  numerical 
solution  of  Madariaga  (1976).  His  calculations  revealed  that 
termination  of  slip  actually  occurs  only  when  healing  phases 
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propagate  inward  from  the  fault  edge  at  the  P  and  S  velocities  of 
the  medium.  The  advantage  of  the  S  &  H  model  Is  that  its 
geometrical  far-field  radiation  can  be  obtained  analytically  in 
closed  form.  Given  the  geometry  shown  in  Figure  36,  the  geometrical 
far-field  acceleration  can  be  written 


(8/a)' 


sin  2e  cos  i  I 


*Ue  =  KTs r  cos  2e  cos  4  h 


(7.1) 


U«f  =  Tnr  cos  8  sin  4  r6 


with 


4*6/  L 

I  =  - H ( t )  -  H ( t  -  y  (1  -  K)) 

c  (1  -  K L)L  V 


where 

6  3  slip  velocity  at  center  of  crack 

o 

V  =  rupture  velocity 

L  =  fault  radius 

K  =  V/C  sin  e 
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t  =  reduced  time 

a,  6  =  compressional  and  shear  velocities,  repectively  in  the 
source  region 

C  =  a  or  8 

H(t)  =  Heaviside  step  function 

ait)  =  delta  function 

The  time  domain  behavior  is  shown  in  Figure  37.  Acceleration 
initiates  as  a  step  which  continues  until  the  arrival  of  the 
stopping  phase  from  the  edge  of  the  fault  nearest  the  observer. 
From  that  time  until  the  time  of  arrival  of  the  stopping  phase  from 
the  farthest  edge  of  the  fault,  the  acceleration  is  constant. 

There  are  several  noteworthy  characteristics  of  the 
accelerations.  The  high  frequency  content  is  clearly  dominated  by 
the  stopping  phases,  which  are  (mathematically)  delta  functions. 
These  cause  the  far-field  displacement  spectrum  to  decay  as  u  at 
high  frequencies.  The  amplitudes  of  the  first  stopping  phase  and 
the  step  initiating  the  motion  contain  terms  of  K(1  -  K)'1  and 
(1  -  K^)  respectively.  These  terms  make  the  amplitudes  of  the 
early  parts  of  the  motion  rather  strong  functions  of  the  rupture 
velocity  and  the  azimuth  of  the  observer  from  the  fault  normal.  For 
example,  with  a  rupture  velocity  of  0.9b,  where  s  is  the  shear  speed 
of  the  medium,  the  amplitudes  of  the  initiation  phase  in  the  plane 
of  the  fault  and  at  the  fault  normal  differ  by  more  than  a  factor  of 
27.  For  a  given  slip  velocity,  changing  the  rupture  velocity  from 
0.8b  to  0.9b  causes  the  initiation  phase  to  increase  by  a  factor  4.5 
in  the  plane  of  the  fault,  but  only  by  1.1  at  the  fault  normal. 

The  delta-function  dependence  in  the  stopping  phase  will,  of 
course,  be  smoothed  by  an  attenuating  medium,  but  even  with 
reasonable  values  for  the  intrinsic  attentuation  Q"*,  it  appears 
that  the  predicted  stopping  phases  are  much  too  strong  to  be 
consistent  with  near-field  observations.  The  amplitude  of  the  phase 
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Figure  37.  Typical  shape  of  the  Sato  and  Hirasawa  model  geomet¬ 
rical  far-field  acceleration.  Arrows  represent 
delta  functions. 


can  be  roughly  estimated,  given  linear  attentuation,  using  the 
asymptotic  formulas  of  Kjartansson  (1979).  He  suggests  the  peak 
time  domain  amplitude  of  a  causal,  attenuated  pulse  to  be  roughly 
CQ/R  times  the  strength  of  the  input  delta  function,  where  C  is  the 
signal  velocity  and  R  is  the  distance  traveled.  For  example,  Hanks 
(1974)  suggested  that  the  initiation  event  for  the  1971  San 
Fernando,  California,  earthquake  had  a  stress  drop  of  350  to  1400 
bars  over  a  fault  radius  of  3  to  6  kilometers.  We  can  estimate  the 
relative  peak  amplitudes  of  the  initiation  and  stopping  phases 
observed  at  Pacoima  dam,  R  =  15  kilometers  and  very  near  the  plane 
of  the  fault.  Figure  38  shows  the  values  predicted  by  the  S  &  H 
model  for  a  400  bar  stress  drop  and  5  kilometer  fault  radius,  for  a 
0  of  100.  Except  when  the  rupture  velocity  is  near  the  shear 
velocity,  the  stopping  phase  is  estimated  to  be  considerably  larger 
than  the  initiation  phase.  Only  relatively  slow  rupture  velocities 
can  predict  reasonable  values.  Since  the  absolute  values  depend  on 
many  poorly  constrained  parameters,  the  relative  values  are  the  most 
meaningful.  If  actual  recordings  are  examined  (Figure  39),  the 
step-like  feature  is  evident  on  all  three  components  of  motion  about 
2.5  seconds  into  the  record.  After  approximately  0.6  seconds,  there 
is  a  sudden  downward  phase  indicating  some  kind  of  stopping,  but  the 
peak  amplitudes  are,  at  best,  comparable  to  the  initiation  phase  on 

the  order  of  0.5  g.  A  better  evaluation  of  the  model  is  achieved  by 

using  the  stopping  time  (t  =  0.6  sec)  to  constrain  the  fault  radius 
L  =  Vt/(1  -  K).  When  this  is  done,  the  predicted  stopping  phase  is 
7  to  8  times  larger  than  the  predicted  initiation  phase  for  all 
rupture  velocities  between  0.5b  and  0.99b. 

Clearly,  the  stopping  phases  of  the  S  &  H  model  are  much  too 
strong,  in  a  linearly  attenuating  medium,  to  be  reasonable.  There 
are  numerous  explanations  for  this  inadequacy  of  the  model,  and  a 
detailed  discussion  will  not  be  undertaken.  There  are  two  items 
worth  discussing  --  the  mode  of  healing  and  the  abruptness  with 
which  rupture  growth  terminates.  The  S  8  H  solution  does  not  treat 

healing  rigorously.  Numerical  solutions  suggest  that  the  healing  of 

the  fault  when  rupture  propagation  terminates  is  not  instantaneous, 
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radius,  assuming  a  Q  of  100. 
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Processed  ground  acceleration  observed  at  Pacoima  Dam  (from  EERL,  1973) 


as  assumed  by  the  S  &  H  model,  but  propagates  inward  at  the  seismic 

signal  velocities.  The  approximate  treatment  of  healing  is  probably 

not  as  influential  in  controlling  peak  acceleration  as  it  may  at 

first  seem;  Madariaga's  (1976)  numerical  solution,  which  treats 

_2 

healing  more  rigorously,  has  an  u  far-field  displacement 

spectrum  indicating  that  it,  too,  gives  singular  acceleration. 

Of  greater  importance  is  the  manner  in  which  rupture  growth 
stops.  In  the  S  &  H  model,  rupture  stopping  is  very  idealized.  The 
propagating  rupture  decelerates  instantaneously  along  a  smooth, 
prescribed  boundary  (a  circle).  The  values  of  observed  peak 
accelerations  suggest  that  this  approximation  is  unacceptable  for 
predicting  the  characteristics  of  high  frequency  radiation. 

Clearly,  a  decelerating  model  of  rupture  is  needed.  The 
D-model  of  Boatwright  (1980)  would  be  an  appropriate  choice,  but  at 
present,  no  simple  far-field  representation  of  the  motion  from  this 
model  is  available.  Although  an  analytic  far-field  solution  is  not 
essential,  it  does  add  considerable  flexibility  to  the  calculations 
and  facilitates  interpretation  of  computed  waveforms.  As  an 
alternative  to  employing  a  rigorous  numerical  solution  for  a 
decelerating  rupture  model,  we  choose,  for  simplicity,  to  alter  the 
form  of  the  far-field  radiation  from  the  S  &  H  model  to  have  a 
deemphasized  stopping  phase.  The  most  convenient  way  to  do  this  is 
to  assume  that  the  rupture  deceleraties  over  an  annular  region  of 
width  S.  The  stopping  phase  will  then  be  smoothed  out  over  a  time 
related  to  S. 

Boatwright's  (1980)  0-model  provides  a  convenient  parameter¬ 
ization  of  the  effect.  He  introduces  a  parameter,  which  we  will 
call  y  (Boatwright's  v  >  1),  which  controls  the  time  of 

the  stopping  process.  For  a  radius  L  and  rupture  velocity  V,  the 
rupture  remains  uniform  out  to  a  radius  yL,  at  which  point  the 
rupture  decelerates.  The  duration  of  the  uniform  rupture  process  is 
then  -Tj-,  and  it  is  assumed  that  the  rupture  reaches  the  boundary 
at  a  time  For  y  =  1.  this  model  is  identical  to  the  S  &  H 

model.  For  0  <  y  <  1,  the  stopping  radiation  will  be  smoothed. 
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One  could  solve  directly  for  the  radiation  due  to  Boatwright's 
0 -model ,  but  since  the  details  of  the  stopping  process  of  such  a 
model  is  arbitrary  anyway,  it  seems  more  convenient  to  instead 
assume  some  simple  form  for  the  characteristics  of  the  stopping 

phases  which  satisfies  all  gross  conservation  conditions  and 
qualitatively  approximates  the  actual  far-field  solution.  The 
far-field  radiation  from  a  dislocation  type  source  must  conserve  the 

integrals  of  acceleration  and  velocity  and  have  an  intergral  of 

displacement  proportional  to  seismic  moment.  As  an  alternative  to 
the  0 -model ,  we  simply  replace  the  delta-function  stopping  phases 
with  steps  of  duration  consistent  with  the  time  of  deceleration  of 
the  D-model  as  perceived  at  the  observer,  and  choose  amplitudes 
which  conserve  the  integrals  of  acceleration,  velocity  and 
displacement.  An  example  is  shown  in  Figure  40. 

In  summary,  the  Sato  and  Hirasawa  crack  approximation  appears 
to  be  unreasonable  when  high  frequency  radiation  is  of  interest. 

The  large  stopping  phases  predicted  by  the  model  are,  in  general, 
too  large  to  be  consistent  with  observed  values  of  acceleration 
close  to  earthquakes.  A  more  acceptable  model  for  the  radiation 
from  a  stress  concentration  apparently  requires  a  deemphasis  of  the 
high  frequency  radiation  due  to  the  stopping  of  rupture.  Here  we 
choose  an  alternative  representation  of  the  far-field  radiation  from 
an  isolated  release  of  stress  which  includes  a  smoothing  of  these 
stopping  phases. 

The  source  parameters  for  the  three  discrete  sources  used  here 
are  summarized  in  Table  11.  The  orientation  angles  are  defined  as 
in  Aki  and  Richards  (1980),  Chapter  4.  The  spatial  coordinate 
system  used  is  centered  at  the  epicenter  of  the  first  source  and 
positive  X  and  Y  is  in  the  north  and  east  directions  respectively. 
The  moment-magnitude  was  determined  using  the  relationships  proposed 
by  Hanks  and  Kanamori  (1979). 

7.4  FREE  SURFACE  CORRECTION 

For  most  near-field  applications,  the  interaction  of  the 
radiation  with  the  free  surface  and  surficial  layering  can  be  quite 
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Figure  40.  Comparison  of  typical  shapes  of  far-field  acceleration  from  the 
Sato  and  Hirasawa  model  and  the  modified  model  used  in  this 
study. 
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Table  11 

SOURCE  PARAMETERS 


#1 

n 

#3 

Fault  strike  (deg) 

-75.0 

-75.0 

-75.0 

Dip  angle  (deg) 

53.0 

29.0 

29.0 

Slip  angle  (deg) 

76.0 

90.0 

90.0 

X-coordinate  (km) 

0.0 

-9.0 

-15.5 

Y-coordinate  (km) 

0.0 

3.0 

4.0 

Depth  (km) 

13.0 

6.0 

2.5 

Fault  radius  (km) 

Fraction  of  radius  for  which 

6.0 

3.0 

1.0 

rupture  is  uniform 

0.4 

0.8 

0.8 

Rupture  velocity  (km/sec) 

2.31 

2.31 

2.31 

Slip  velocity  (cm/sec) 

600.0 

200.0 

200.0 

Initiation  time  (sec) 

0.0 

7.4 

11.6 

Moment  magnitude 

6.4 

5.6 

4.6 
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important  to  determining  the  characteristics  of  the  surface  ground 
motion.  For  the  case  of  the  San  Fernando  earthquake,  the  medium  in 
the  vicinity  of  the  event  can  be  adequately  approximated  by  an 
elastic  half-space,  and  surficial  layering  does  not  appear  to  be 
important  in  the  epicentral  region.  There  are  a  number  of  exact 
algorithms  for  determining  the  free  surface  motions  from  point 

dislocation  sources  (Johnson,  1974,  for  example).  These  algorithms 
are  not  expensive  when  the  response  for  few  source-receiver  pairs 
are  of  interest.  Standard  methods  for  computing  the  atmospheric 
coupling  must  evaluate  a  surface  integral  of  the  ground  acceleration 
for  each  time  desired  in  the  acoustic  signal.  For  this  case,  use  of 
rigorous  free  surface  responses  would  be  much  too  costly. 

As  an  alternative  to  complete  free  surface  calculations,  we 
employ  an  approximate  method  first  suggested  by  Knopoff  et  al . , 
(1957).  The  total  response  is  approximated  by  the  response  of  an 
incident  plane  wave  with  incident  angle  corresponding  to  the 

geometrical  arrival  in  the  complete  response.  The  correction  factor 
is  a  complex  scale  factor  which  can  be  obtained  analytically.  The 
scale  factor  is  applied  such  that  the  output  signal  is  the  input 
signal  times  the  real  part  of  scale  factor  plus  the  H  lbert 
transform  of  the  input  signal  times  the  imaginary  part  of  the  scale 
factor.  The  final  response  will  not  include  Rayleigh  waves  or  S  to 
P  phases  refracted  along  the  surface. 

An  accurate  derivation  of  the  method  can  be  found  in  Anderson 
(1976)  (Knopoff  et  al . ,  1957,  contains  some  misprints).  Given  a 
P-wave  with  incident  angle  i  or  an  S-wave  with  incident  angle  j,  we 
define  a  ray  parameter: 

.  _  sin  i  _  sin  j 

5  "  a  "  B 

For  incident  P-waves,  the  scale  factor  is  always  real  and  can  be 
expressed  as: 
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where 


2 

n  -  2s  cos  i  sin  2j  +  cos  2j 

a  ~7~  • 

rP  (RP)  is  the  ratio  of  the  predicted  vertical  (radial) 
response  to  the  whole  space  vertical  (radial)  response. 

For  incident  SH  waves,  the  transverse  scale  factor  is  two, 
independent  of  incident  angle.  For  SV  waves,  the  scale  factor  is 
real  for  sub-critcal  incidence  (sin  j  <  s/a)  and  complex  for 
super-critical  incidence.  For  subcritical  case. 


rs  =  4s  cos  i  cot  j 
v  ciD  *" 

0s  _  2  cos  2j 

H  '  “77^  * 

For  super-critical  incidence, 


4  sb  cot  j 

0, 


( i  ( p  -  J) 


R 


s 

H 


2  cos  2j 
62D‘ 


( i  ( p  +  e ) 
e 


where 


b  =  ,  /$2  -  1/a2 

V 


,,  ,  4 


)  +  (2sb  sin  2j) 


1/2 


tanp  =  gg.  _sbjinJl 
cos^  2j 


and 


© 


0  sin-1  ^<JiJ 

ir  .  .  .  it 

ir  *  <  J  <  1 


The  approximation  used  here  can  be  summarized  as  a 
high-frequency,  first-motion  approximation.  It  is  exact  for  SH 
waves  and  excellent  for  incident  P-waves  at  all  incident  angles  and 
for  subcritical  S V  waves.  Near  the  critical  angle,  it  is  poor. 
Beyond  critical,  it  performs  quite  well  at  times  near  the  arrival 
time  of  the  geometrical  arrival,  but  fails  at  times  where  a  Rayleigh 
wave  or  an  S-P  converted  phase  is  expected.  Anderson  (1976)  shows 
several  comparisons  between  the  response  predicted  by  the  method 
given  here  and  complete  half-space  and  whole  space  responses. 


7.5  DESCRIPTION  OF  ROUTINES 

One  master  routine,  SFVERT,  and  five  slave  routines  generate 
approximate  vertical  accelerations  for  a  simulation  of  the  1971  San 
Fernando,  California  earthquake.  The  version  summarized  here 
returns  the  acceleration  for  a  given  location  at  a  given  time.  It 
does  not  generate  a  complete  time  history  at  each  call.  The 
routines  are  accessed  by: 
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CALL  SFVERT  (DT,  XR,  YR,  XC,  T) 


where 

OT  =  Input  time  increment  used  only  for  smoothing  the 
response  to  prevent  aliasing.  The  value  used  should 
be  comparable  to  the  time  increment  of  a  final 
seismogram  or  pressure  signature  desired. 

XR,YR  =  Input  cartesian  coordinates  of  desired  receiver  in 
kilometers.  The  coordinate  system  is  positive  to 
the  north  and  east  and  centered  at  the  earthquake 
epicenter. 

2 

XC  ■  Ouptut  vertical  acceleration  in  cm/sec  ,  positive 
down. 

T  =  Input  time  of  desired  response. 

Each  call  of  SFVERT  initiates  the  following  sequence  of 

operations: 

1.  The  vector  pointing  from  the  source  to  the  receiver  is 

rotated  from  the  global  coordinate  system  to  a  source 

system. 

2.  Radiation  pattern  is  computed  and  a  motion  vector  and 

far-field  time  history  is  obtained. 

3.  The  motion  vector  is  rotated  back  to  the  global 

coordinate  system. 

4.  The  free  surface  correction  is  computed. 

5.  The  result  is  convolved  with  the  far-field  time  history. 

Steps  (2)  through  (5)  are  repeated  for  P,  SV  and  SH  motion.  The 

total  source  model  consists  of  three  individual  sources,  and  all 

steps  are  repeated  for  each  individual  source. 

The  routines  performing  the  tasks  are: 

SFVERT  -  Master  routine,  contains  source  data,  generates 

far-field  time  histories,  applies  free  surface 
correction. 
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SFVERI 


Performs  coordinate  transformations,  computes 
radiation  patterns. 

FSURF  -  Computes  P  and  SV  free  surface  correction. 

ROTATE  -  Performs  vector  rotations. 

FSMOO  -  Smoothing  routine. 

DETER  -  Assists  in  computation  of  time  histories. 

The  time  histories,  before  the  free  surface  correction  is 
applied,  consist  of  a  series  of  steps  in  acceleration.  They  will, 
in  theory,  have  zero  area  and  double  area  proportional  to  seismic 
moment.  In  practice,  these  conditions  are  not  met  exactly  due  to 
discretization.  The  noncasual  free  surface  correction  also  affects 
the  areas  in  practice. 

The  source  parameters  for  each  crack  are  specified  in  a  DATA 


statement  in 

the 

master  routine  SFVERT.  The  individual  parameters 

are  as  follows: 

STRI 

Fault  strike  (degrees). 

DIPS 

= 

Fault  dip  (degrees). 

SLIPS 

= 

Slip  angle  (degrees). 

XS,YS 

= 

Cartesian  coordinates  of  the  crack  epicenter  (km) 

HS 

= 

Source  depth  (km). 

LS 

a 

Crack  radius  (km). 

GAMS 

a 

1  1/2 

See  Boatwright's  paramter  (— — ) 

DOS 

a 

Slip  velocity  at  crack  center  (cm/sec). 

VS 

a 

Rupture  velocity  (km/sec). 

TS 

= 

Origin  time  (sec). 

Orientation  angles  defined  in  the  usual  way.  See  Aki  and  Richards 
( 1980) ,  Chapter  4. 
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7.6  COMPARISON  WITH  NEAR-FIELD  OBSERVATIONS 

In  this  section  we  compare  the  near-field  response  predicted 
by  the  procedures  described  above  with  the  observed  ground  motion  at 
Pacoima  Dam.  The  accelerometer  at  Pacoima  Dam  is  one  of  over  90 
strong  motion  instruments  which  triggered  during  the  San  Fernando 
earthquake.  It  was  located  only  a  few  kilometers  from  the  rupture 
surface  and  at  an  epicentral  distance  of  about  8  to  14  kilometers, 
depending  on  the  reference  cited  (Figure  41).  The  next  nearest 
recording  of  the  ground  motion  was  located  more  than  20  kilometers 
from  the  epicenter.  The  motion  recorded  at  Pacoima  Dam  has  been 
quite  controversial .  When  recorded  in  1971,  the  horizontal 
accelerations  recorded  were  the  largest  known  to  date.  The 
instrument  was  located  on  the  side  of  a  ridge,  and  it  was  one  of  the 
first  instances  where  topographic  effects  were  considered  as 
important  to  the  motion  recorded.  After  the  event,  the  instrument 
was  found  to  be  tilted  from  its  preearthquake  position  (Trifunac  and 
Hudson,  1973).  Though  debate  still  occurs  as  to  the  credibility  of 
the  high  frequencies  recorded  (empirical  studies  of  ground 
acceleration  often  ignore  the  recording),  the  intermediate  frequency 
content  appears  reasonable. 

For  simulation  of  acoustic  radiation,  the  vertical  ground 
acceleration  is  needed.  The  resulting  pressure  signatures  in  the 
linear  part  of  the  atmosphere  most  resemble  vertical  ground  velocity 
in  frequency  content.  Though  the  algorithm  described  in  previous 
sections  does  generate  ground  acceleration,  the  acceleration 
produced  is  intended  to  be  reasonable  only  in  the  intermediate 
frequency  range  (.5  to  3  Hz)  which  dominates  ground  velocity  (and 
atmospheric  pressures).  For  this  reason  the  fairest  test  of  our 
methods  is  a  comparison  of  the  recorded  ground  velocities  at  Pacoima 
Dam,  and  velocity  obtained  through  integration  of  the  output 
predicted  by  our  algorithm. 

Figure  42  shows  the  ground  velocity  at  Pacoima  Dam  with  that 
simulated  with  our  procedures.  Given  the  simplicity  of  our  source 
and  propagation  models,  the  comparison  is  quite  favorable,  and  most 
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Source  =2 


Pacoima 
Dam  A 

H*  Source  -3 


Location  of  Pacoima  Dam  accelerograph  site  with  the  epicenters 
of  the  three  discrete  sources. 


Comparison  of  observed  and  computed  vertical  ground  velocities  at 
Pacoima  Dam. 


previous  attempts  to  model  the  ground  velocity  are  not  much  more 
effective.  The  major  deficiency  concerns  the  response  due  to  the 
second  source  of  radiation.  Our  simulation  contains  approximately 
the  proper  amount  of  energy,  but  appears  to  be  phase  shifted  by 
about  90  degrees.  Past  modeling  efforts  have  not  had  success  at 
dealing  with  this  arrival  either  (see  Boore  and  Zoback,  1974,  and 
Bache  and  Barker,  1978).  Our  methods  have  some  restrictions  which 
may  be  responsible  for  the  discrepancy.  Other  models  of  the  source 
of  radiation  located  above  the  kink  of  the  fault  have  suggested  that 
the  rupture  may  have  been  unilateral.  Our  source  model  is  limited 
to  bidirectional  rupture.  The  incident  geometrical  arrivals  at 
Pacoima  Oam  from  the  source  in  question  is  very  near  the  critical 
angle  for  incident  S-waves.  The  poor  performance  of  the  first 
surface  correction  used  here  near  critical  incident  angles  may  also 
contribute  to  the  poor  fit. 

In  summary,  we  have  presented  a  simple  algorithm  for  the 
simulation  of  the  near-field  ground  accelerations  from  the  1971  San 
Fernando  earthquake  which  shows  reasonable  agreement  with  observed 
ground  motion  in  tire  intermediate  frequency  band.  Our  intent  was 
not  to  develop  a  better  model  of  the  earthquake,  but  to  use  the 
results  of  previous  studies  and  to  apply  a  few  carefully  chosen 
approximations  to  develop  an  inexpensive  tool  for  providing 
approximate  near-field  ground  accelerations  which  can  be  used  for 
the  simulation  of  the  atmospheric  disturbances  expected  for  a 
moderate  sized  earthquake. 
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APPENDIX  A.  RUPTURE  MODEL  FOR  MULTI-MATERIALS  VERSION  OF  I4TRES 


This  appendix  describes  the  algorithm  modifications  required 
to  incorporate  the  twc-degree-of-freedom,  slip-weakening  rupture 
model  (Day,  1979)  into  the  14TRES  code.  The  notation  basically 
follows  that  of  the  Advanced  Computing  Laboratory  document  ACL-80-035. 
Modifications  are  confined  to  the  fault  plane  equations. 


On  the  slip  plane,  we  permit  motion  in  the  u  component,  and 
form  u  and  u'  bv  analogy  with  v  and  v': 

u(JN00£,KN00E,L$LIP)  » 

l/2[u+( JNODE, KNODE, LSI IP) +u""( JNODE, KNQDE.LSL IP )] 
u 1 (JNODE, KNOOE.LSLIP )  * 

l/2[u+( JNODE, KNODE.LSL IP)  -u~( JNODE.KNODE.ISLIP) ] , 

with  similar  definitions  for  Q,  jj,  u ' ,  and  u'.  Likewise,  we  form  'HGACCX  and 
HGACCX 1  by  analogy  with  HGACC  and  HGACC': 

For  JNODE  *  JSLMN..JSLMX  ( JSEGMN<JSLMN,JSEGMX>JSLMX) 

(HGACCX  »  )  HGACC ( JNODE, 1)  -  4x(JN0DE)*4y(KN0DE)*.5* 

*| MINV( JNODE-1/2 , KNODE-1/2 .IN00E-1/2 )* 

*[ 62 ( LNODE-1/2 ) *TQX2 ( JNODE-1/2 .KNODE-l /2 ,L  NODE-1 /2 ) +«z (LNODE-1/2 ) 
*TQX3( JNODE-1/2 .KNODE-1/2 , LNODE-1/2 ) ] 

+MINV( JN00E+1/2 .KNODE-1/2 .LNODE-1/2 }* 

*[-«  (LNODE-1 /2 ) *TqX2 ( JN00E+1 /2 , KNODE-1/2 , LNODE-1/2 )+6Z (LNODE-1 /2 ) 
*TQX3( JNOOE+1/2 , KNODE-1/2 , LNODE-1/2 )] 

+MINV( JNQ0E-1/2 ,KN0DE+l/2 , LNODE-1/2)* 

*[«z( LNODE-1 /2)*TQX2( JNODE-1/2, KNODE+1/2, LNODE-1/2 )-«z (LNODE-1/2) 
*TQX3( JNODE-1/2 ,KN0DE+l/2 , LNODE-1 /2 )] 

+MINV( JN00E+1/2 ,kN0DE+l/2 .LNODE-1/2 )* 

*[ -62 (LNODE-1 /2 )*TQX2 ( JN0DE+1/2 .KNODE+l / 2 , LNODE-1 / 2 ) -5Z ( LNODE-1/2 ) 
*TQX3( JN00E+1/2 .KNODE+1/2 , LNODE-1/2 )] 
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+MINV( JNODE-1/2 , KNODE-1/2 , LNODE+1 /2 )*[2*?2 (LNODE )*TQX1 (JNQDE-1/2 ,  KNODE-1/2) 
-az ( LNODE+1 / 2 ) *TQX2 ( JNOOE-1 / 2 , KNODE-1 12 , LNODE+1/2 ) — «2 ( LNODE+1 12 ) 

*TQX3( JNODE-1/2, KNOOE-1/2, LNODE+1/2)] 

+MINV(JN0DE+l/2, KNOOE-1/2, LN0DE+l/2)*[-2*;2(LN0DE)*TQXl( JNODE+1/2, 

KNOOE-1/2) 

+«z ( LNODE+1 /2 )*TQX2( JNODE+1/2 ,KNODE-l/2 ,LN0DE+l/2 )-dz(LNQDE+l/2 ) 

*TQX3( JNODE+1 /2 ,KNOD£-l /2 , LNODE+1 /2 )] 

+MINV( JNODE-1/2 ,KN00E+l/2 ,LNOOE+l/2 )*[-2*Sz(LNQOE)*TQXl (JNODE-1/2 , 

KNODE+1/2) 

-«z(LN00E+l/2)*TQX2( JNODE-1/2, KNODE+1/2, LNODE+1 /2 ) +6Z ( LNODE+1 /2 ) 

*TQX3 ( JNODE-1 /2 , KNODE+1/2 , LNODE+1 /2 ) ] 

+MINV(JN00E+l/2,KN00£+l/2,LN00E+l/2)*[2*S2(LN0DE)*TQXl(JN0DE+l/2,  KNODE+1/2) 
+«z (LNODE+1 12 )*TQX2 (JNOOE+1 /2 ,KNODE+l /2 , LNODE+1 /2 ) +«z (LNODE+1 / 2 ) 

*TQX3 ( JNOOE+1 /2, KNODE+1/2, LNODE+1/2)]  [. 


Similarly, 

HGACCX' (JNODE)  =  5 x(JNODE)*?y(KNODE)  *  .5* 

* | MINV( JNODE-1/2 , KNODE-1 /2,LN00£-i/2 )* 

*[-dz(LN0DE-l/2 )*TQX2 (JNODE-1/2 ,KNODE-l/2 .LNODE-1/2 )-«z (LNODE-1/2 ) 
*TQX3 (JNODE-1 12 ,KN0DE-1 /2 , LNODE-1 / 2 ) ] 

+MINV( JNODE+1/2 .KNODE-1/2 ,LN0DE-l/2 )* 

*[ sz ( LNODE-1 /2 )*TQX2( JNOOE+1 12 , KNODE-1 / 2 , LNODE-l 12 )-6Z ( LNODE-1 /2 ) 
*TQX3(JN00E-l/2, KNODE-1/2, LNODE-1/2)] 

+MINV( JNOOE-1/2 ,KNODE+l/2 , LNODE-1 /2 )* 

*C-«z (LNODE-1/2 )*TQX2(JNODE-l/2, KNODE+1/2, LNODE-1 /2)+6Z (LNODE-1/2) 
*TQX3 (JNODE-1 12 .KNODE+l 12 , LNODE-1 12 ) ] 

+MINV( JNODE+1/2 , KNODE+1/2 ,LN00E-l/2 )* 

*[ 6Z (LNODE-1 /2 )*TQX2 (JNODE+1/2 .KNODE+l 12 , LNODE-1/2 ) +az ( LNODE-1 12  ) 
*TQX3( JNODE+1/2 , KNODE+1/2 , LNODE-1/2 ) ] 

+MINV( JNODE-1/2, KNODE-1/2 .LNODE+1/2 )*[2*^Z(LN0DE)*TQXID( JNODE-1/2 , 
KNODE-1/2) 

-5Z(LN0DE+1/2)*TQX2 (JNODE-1/2, KNODE-1/2, LNODE+1/2 )-«z (LNODE+1/2) 
*TQX3( JNODE-1/2, KNOOE-1/2, LNODE+1/2)] 
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* 


+MINV(JN0DE+l/2, KNODE-1/2, LNODE+l/2)*[-2*5z(LNODE)*TQXID(JNODE+l/2, 

KNODE-1/2) 

+az(LN00E+l/2)*TQX2(JN0D£+l/2, KNODE-1/2, LN0DE+l/2)-az(LN0DE+l/2) 

+TQX3 ( JNQDE+1/2 ,KN00E-l/2 .LNODE+1/2 ) ] 

+MINV( JNOOE-1/2 ,KNODE+l/2J.NODE+l/2 )*[-2*Sz(LN0DE)*TQXID( JNODE-1 /2 , 

KNODE+l /2) 

-6z(LN00E+l/2 )*TQX2( JNOOE-1/2 ,KN00E+l/2 ,LN0DE+l/2  )+<Sz(  LNODE+1/2 ) 

*TQX3 ( JNODE-1 /2 , KNODE+l 12 , LNODE+1 / 2 ) ] 

+MINV(JNODE+l/2, KNOOE+1/2, LNODE+l/2)*[2*;z(LNQDE)*TQXID(JNODE+l/2, 

KNOOE+1/2) 

+«z ( LNODE+1 12 ) *TQX2 ( JNOOE+1 /2 , KNODE+l / 2 , LNODE+1 12 )+6Z ( LNOOE+1 12 ) 

*TQX3( JNOOE+1 /2, KNOOE+1/2, LNODE+1 /2)] j . 

In  the  above,  TQXl  and  TQXID  are  defined  by 

TQX1(JCELL, KCELL)  »  p(JCELL, KCELL, LCELL)*BR(JCELL, KCELL, LCELL) 

*6z{ LCELL ) / [s  z ( LCELL+1/2 ) *«x ( KCELL ) ]*HG1 ( JCELL , 1 ) 

TQXID(JCELL.KCELL)  -  p(JCELL, KCELL, LCELL)*BR(OCELL, KCELL, LCELL) 

*6Z  ( LCELL )  /[c;z(  LCELL-1/2  )*6x  (KCELL ) ]*HG1X ' (JCELL ) , 

with 

HG1( JCELL, 1)  *  —Q( JCELL— 1 / 2 ,KCELL-l/2 , LCELL-1/2 ) +iJ( JCELL+1 /2 , KCELL— 1 / 2 , LCELL— 1 / 2 ) 
+Q( JCELL-1 12 .KCELL+l /2 ,  LCELL-1/2  )-IT(  JCELL+1/2 .KCELL+1/2 , LCELL-1/2 ) 
HG1X 1 ( JCELL)=  -Ci  ‘ (OCELL-1/2 .KCELL-1/2 .LCELL-l /2 )+u ‘ (JCELL+1/2, KCELL-1/2, LCELL-1/2) 
+u'( JCELL-l/2,KCELL+l/2,LCELL-l/2)-u( JCELL+1/2, KCELL+1/2, LCELL-1/2). 

Off  the  slip  plane,  TQX1D(JCELL, KCELL)  *  0. 

Next  we  form  the  acceleration  terms  XSTART  and  YSTART: 

XSTART(JNODE)  -  Q.5*«z(LN0DE-l/2)*Sy(KN0DE )*j -MXINV( JNODE, KNODE-1/2 .LNODE-1/2 ) 

*C® w( JNOOE+1/ 2, KNODE-1/2, LNODE-1/2 )-o( JNODE-1 /2, KNODE-1/2 ,LN0DE-l/2 )] 

-MX INV (JNODE , KNODE+l /2 .LNOOE-l / 2 ) 
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*[a  (JNODE+1/2 , KNODE+l /2 , LNODE-1/2 )-o  (JNODE-1/2* KNOO£+1/2,L NODE-1 /2)] | 

+.5*«z( LNODE+1/2 , )*Sy(KNOOE)*j  MXINV( JNODE,KNOOE-l/2,LNQOE+l/2) 

*[o  (JNQOE+1/2 , KNODE-1/2 , LNODE+1/2 )-<*Yy(  JNODE-1/2,  KNODE-1/2,  LNODE+1/2 )] 
+MX INV ( JNODE , KNODE+l / 2 , LNOOE+1 / 2 ) 

*[a  v (JN0DE+l/2,KN0DE+l/2,LN0DE+l/2) -a  (JNODE-1/2, KNODE+l /2, LNODE+l /2)][ 
+. 5*«z (LNODE-1 /2 ) *£x (JNODE ) *  j -MY INV ( JNOOE-1 /2 ,LNODE-l/2) 

*[axy( JNODE-1/2, KNODE+1/2, LNODE-1/2 )-oxy( JNODE-1/2, KNOOE-1/2, LNODE-1/2)] 
-MY INV ( JNODE+l /2 , LNODE-1 /2 ) 

*|>xy( JNODE+1/2 ) ,KN0D£+l/2 , LNODE-1/2 )-axy( JNODE+1/2 , KNODE-1/2 , LNODE-1/2 )] 
+.5*6Z(LN0DE+1/2)*CX( JNODE )*|mYINV( JNODE-1/2 , LNODE+1/2 ) 

*l>xy(JNODE-l/2, KNODE+1/2, LNODE+1/2 )-axy( JNODE-1/2, KNODE-1/2, LNODE+1/2)] 
+MY IN V ( JNOOE+1 /2 , LNOOE+1 / 2 ) 

*£°y(  JNODE+1/2, KNOOE+1/2, LNOOE+1 /2)-axy( JNODE+1/2, KNODE-1/2, LNODE+1/2)]  | 

LNOOE+1 /2  KNODE+l /2  JNOOE+1 /2 

+.5*Sx(JNODE)*Cv(KNODE)*I  I  I 

J  L-LNOOE-1/2  K.KNODE-1/2  J-JNODE-1/2 

[MINV(J,K,L)+(  cxz(J,K,L)  +  yxz  (JNODE))  ] 

+HGACCX'( JNODE). 

YSTART(JNODE)  =  . 5*62 ( LNODE-1 / 2 ) *^y ( KNODE ) * j -MX  I NV ( JNODE, KNODE-1/2 , LNODE-1/2 ) 
*[oxy( JNODE+1/2 .KNOOE-1/2 , LNODE-1/2 )-axy( JNODE-1/2 .KNODE-1/2 , LNODE-1 IZ ) ] 
-MX INV ( JNODE , KNODE+l /2 , LNQDE-1 / 2 ) 

*1>  xy( JNODE+1/2, KNODE+1/2, LNOOE-l/2)-oxy( JNODE-1/2, KNODE+1/2, LNODE-1/2 )]j 
+. 5*5z( LNODE+l /2)*?y(KN00E)*} MX INV( JNODE, KNODE-1/2 ,LNODE+l/2) 
*I>xy(JNODE+l/2, KNODE-1/2, LNQDE+l/2)-axy( JNODE-1/2, KNODE-1/2, LNODE+1/2)] 
+MX INV( JNODE , KNODE+l /2.LNQDE+1/2 ) 

*[a  (JNODE+1/2, KNOOE+1/2 .LNODE+1/2 )-a  ( JNODE-1/2, KNODE+l /2.LNODE+1/2)]  j 

+.5*«z(LNOOE-l/2)*Cx(JNODE)*j-MYINV( JNODE-1/2, LNODE-1/2) 

*[ayy( JNODE-1/2 ,KN0DE+l/2 .LNOOE-1/2 )-oyy ( JNODE-1 /2 .KNODE-1/2 .LNODE-l /2 ) ] 
-MYINV ( JNOOE+1 /2, LNODE-1/2) 

*C®yy( JNODE+1/2, KNODE+1/2, LNODE-1/2 )-ayy( JNODE+1/2, KNODE-1/2, LNODE-1/2)]  j 
+  . 5*5z( LNODE+1/2 )*^x( JNODE)*J  MYINV (JNODE-1/2 .LNODE+1/2) 

*[ayy( JNODE-1/2 , KNODE+1/2 .LNODE+l /2 )-oyy ( JNODE-1/2 , KNODE-1/2 .LNODE+l /2 ) ] 
+MYINV( JNOOE+1 /2 , LNODE+1 / 2 ) 
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*[°yy( JNOOE+1/2 .KNOOE+1/2 ,LN0DE+l/2 )-cyy( JNODE+1/2 .KNODE-1/2 .LNOOE+1/2 ) ] [ 

LNOOE+1/2  KNODE+1/2  JNODE+1/2 

+. 5*Cv ( JNOOE )*SV (KNODE )*£  I  Z 

L.LN00E-1/2  KaKNODE-1/2  J-JNODE-1/2 

[MINV(J,K,L)*  (  oyz(J,K,L)  +  Yyz  (JNOOE)  )  ] 

+HGACC’ (JNOOE), 

in  which  y  and  y  are  components  of  shear  prestress. 

Jr  L  A  L 

Finally,  we  form  a  quantity  A  at  each  slip  plane  node,  defined  by 

LNODE+1/2  KNODE+1/2  JNODE+1/2 

A(JNODE)  -  >5*Sx( JNOOE )*Sy(KNODE )*£  2  2  [MINV( J,K,L)] 

L-LN00E-1/2  K-KNODE-1/2  J-JNODE-1/2 

The  following  input  parameters  are  required  for  the  rupture  model:  SFRAC, 

SMAX,  SMIN,  RUPV,  RCRIT. 

Then,  at  each  node  on  the  slip  surface,  at  time  TTIME,  we  determine  the  nodal 
Quantities,  SIGXZS,  SIGYZS,  YZBFR,  YZSTAR,  ISL  (omitting  the  subscripts  "JNODE, 
KNODE" )  according  to  the  following  algorithm: 


D 

STRF 

R 


[(u1)2  +  (v')2]  1/2  if  [  (u '  )2  +  (J')2]  >  0 

(XSTART2  +  YSTART2)1/2  if  [ (u  1  )2  +  (v’)2]  ,  0 

u'/s  if  [ (u '  )2  +  (v')2]  >  0 

XSTART/ s  if  C(u')2  +  (v1)2]  -0 


.  (1.  -  D)*SMAX  +  D+SMIN 
-  distance  of  node  from  focus 
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’«■/ s  if  C (u •  )2  +  (v')2]  >  0 

YSTART/s  if  [(u')2  +  (v1)2]  »  0 


TR  .  R/RUPV 


YZ8FR 
ISL  -  1 


(XSTART2  +  YSTART2)1^/! 


If  (TR  .  GT  .  TTIME  .  OR  .  R  .  GT  .  RCRIT)  go  to  10 
F  -  min[(TTIME-TR)/flO.  *  «t)  ,  1.0] 


STRF  2  -  (1.  -  F )*YZBFR  ♦  F*SMIN 
IF  (STRF2.  LT.  STRF)  STRF  =  STRF2 
ISL  =  2 


10  continue 

B  -  .5*  «t*STRF*A 

CX  -  u '  +  «t*XSTART  -  GX*B 

CY  «  v'  +  «t*YSTART  -  GY*B 

C  *  (CX2  +  CY2)1/2 

CX  -  CX/C 

CY  *  CY/C 

If  (C  .  GT  .  B)  go  to  20 

YZSTAR  .  [(v'/«t  +  YSTART)2  +  (u'/«t  +  XSTART)2]1/2/ A 
ISL  »  -ISL 
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SIGXZS  -  (u'/«t  +  XSTART)/A 


SIGYZS  .  (v'/«t  +  YSTART)/A 
Go  to  220 

20  YZSTAR  -  .5*STRF*[(GX  +  CX)2  +  (GY  +  CY)2]1/2 
SIGXZS  -  ,5*STRF*(GX  +  CX) 

SIGYZS  -  .5*STRF*(GY  +  CY) 

220  continue 

Finally,  we  form  the  mean  and  differential  accelerations.  The 
differential  accelerations  are  determined  by 


u'  -  XSTART  -  A  *  SIGXZS 
v'  =.  YSTART  -  A  *  SIGYZS, 

and  the  mean  accelerations  are  determined  by 

(0-)u(JNOOE,l)«  .5*5z(lN00E-l/2 )*  Cy(KNODE)*|MXINV( JNODE, KN0DE-1/2,LN0DE -1/2) 
*[axx(JN00E+l/2, KN0DE-l/2,LN00E-l/2)-oxx(JN00E-l/2, KNOOE-1/2, LNODE-1/2)] 
+MX INV ( JNODE ,  KNODE+l 12 , INODE-1 / 2 ) 

*[a„v(JN00E+l/2 .KNODE+1/2 ,LN00E-l/2 )-o  ( JNODE-1/2 ,KN0DE+l/2 ,LN0DE-l/2]j 

XX  l  XX  1 

+. 5*SZ (LNODE+1 /2)*Sy(KN0DE)*| MX INV ( JNODE, KN0DE-l/2,LN0DE+l/2) 

*[o  ( JNODE+1/2 .KNOOE-1/2 ,LN00E+l/2 )-«xx( JNODE-1/2 .KN00E-1/2 ,LN0DE+l/2 ) ] 

+MX INV (JNODE ,KN0DE+1 12 ,IN0DE+1 /2 ) 

*[ovv( JNODE+1/2, KN0DE+l/2,LN0DE+l/2)-ov  (JNODE-1/2, KNODE+l /2.LNODE+1/2)]  | 

XX  |  XX  j 

+.5*«y(LNOOE-l/2)*Sx(JNODE)*jMYINV(JNODE-l/2,LNODE-l/2) 

*[<Jxy( JNODE-1/2, KN00E+l/2,LN00E-l/2)-oxy(JN00E-l/2, KNOOE-1/2, LNODE-1/2)] 
+MY I NV ( JNOOE+1 12 , LNODE-1 / 2 ) 

*!>xy( JNODE+1/2 ,KNQDE+l/2 .LNODE-1/2 )-oxy( JNODE+1/2 ,KN0DE-l/2 , LNODE-1/2 ) ] j 
+. 5*4z(LNODE+1/2)*Cx(JNODE)*|mYINV( JNODE-1/2, LNODE+1/2) 
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*[<J¥w(v!N00E-1/2  .KNODE+1/2  ,LNQ0E+l/2  )-a  ( JNOOE-1/2  ,KN00E-l/2  .LNODE+1/2 )] 

ajt  aj 

+MY INV ( JNODE+1 /2 ,  LNOOE+l /2 ) 

*Cct¥1/(ONOOE+1/2  ,KN0DE+l/2  .LNODE+1/2  )-<?  (  JNODE+1 /2  .KNODE-1/2  .LNODE+1/2 )]! 


KNOOE+1/2 

.5Cx(JNODE)*Cy(KNODE)*2: 

K=KN0DE-l/2 


JNODE+1 /2 

Z 

JaJNODE-1/2 


[MINV( J.K, LNOOE+l /2 )*<j  (  J.K,  LNODE+1/2 )-MINV( J.K.LNODE-1/2 )*o  ( J.K, 

LNODE-1/2)] 


+  .5*4x(JNODE)*i;y(KNOO£)*  (  SIGXZS(JNOOE)  -  Yxz  ( JNODE )  )* 


KNODE+1/2 

z 

K=  KNODE— 1 / 2 


JNODE+1 /2 

V 

J=JN0DE-l/2 


[MINV(J,K,LN0DE-l/2)-MINV( J.K, LNODE+1/2)] 
+HGACCX( JNODE, 1) , 


(V«)v( JNODE, 2)-  . 5 6Z (LNODE-1/2 ) *Cy (KN00E)*| MX  IN V( JNODE, KNODE-l/ 2, LNODE-1/2 ) 
*Oxy(  JN00E+1/2 .KNODE-1/2 .LNODE-1/2 )-oxy ( JNODE-1/2 , KNODE-1/2 .LNODE-l / 2 ) ] 
+MX INV ( JNODE , KNOOE+1 /2 , LNOOE-1 / 2 ) 

*Cay v( JNODE+1 /2 .KNODE+1/2 .LNODE-1/2 )-°Xw( JNODE-1/2 .KNODE+1/2 .LNODE-1/2 )  ] ! 

*7  I  Ay 

+.5*«z(LN0DE+l/2)*G(KN00E)*jMXINV( JNODE, KNODE-1/2, LNODE+1/2) 
*f>xy(JNOOE+l/2 .KNODE-1/2 .LNODE+1/2 )-o  ( JNODE-1/2 , KNODE-l /2 .LNOOE+l /2 ) ] 

+MX INV( JNODE .KNODE+l /2 , LNODE+1/2 ) 

*C«rvx/( JNODE+1 /2 .KNODE+1/2, LNOOE+l /2 )-o  ( JNODE-1/2 .KNODE+1/2.LNODE+1/2 )]  J 

+. 5*«z (LNOOE-1 /2 ) *^x( JNODE )*j MY INV (JNODE-1/2 .LNODE-1/2) 

*[ayy ( JNODE-1 12  .KNODE+1/2 .LNODE-1/2 )-ayy ( JNODE-1/2 , KNODE-1/2 .LNODE-l J2 ) ] 
+MYINV( JNODE+1 /2, LNODE-1/2) 

*[ayy(JNODE+l/2, KNOOE+1/2, LNODE-1/2 )-ayy(JNQDE+l/2, KNODE-1/2, LNODE-1/2)] j 
+.5*dz(LNODE+l/2 )*^x( JNODE )*j  MYINV( JNODE-1/2 , LNODE+1/2 ) 

*[ayy( JNODE-1/2 .KNOOE+1/2 .LNOOE+1/2 )-cyy( JNODE-1 /2 .KN00E-1 / 2 , LNODE+1/2)  ] 
+MY INV ( JNODE+1 12 , LNOOE+l / 2 ) 

*[ayy { JNODE+1 /2 .KNODE+l /2 .LNODE+l /2 )-ayy ( JNODE+1 /2 , KNODE-l /2 , LNODE+1 /2 ) ] J 

KNODE+1/2  JNODE+1 /2 

+. 5*SX ( JNODE )*?y(KNODE*£  Z 

K.KN00E-1/2  J=JNODE-l/2 
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[MINV( J,K,LN0DE+l/2 )*ayz( J,K,LN00E+l/2 )-MINV( J ,K,LNODE 
LNOOE-1/2 )] 

+  .5nx(JN00E)*Cy(KN00E)*  (  SIGYZS(JNODE)  -  Yyz  (jnODE)  ) 

JNODE+1/2 

I 

J-JNOOE-1/2 

[MINV(J,K,LN00E-l/2)-MINV(J,K,LN0DE+l/2)] 
+h5SUc(jnode,2) . 
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l/2)*ayz(J,K, 


KN00E+1/2 

*1 

K-KNQDE-1/2 


MMaii 


i aMflUi 
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APPENDIX  B 

I4TRES  TEST  PROBLEMS 


1.  TEST  PROBLEM  I:  BURIED  DISLOCATION  IN  A  UNIFORM  ELASTIC 
HALFSPACE 

In  this  problem,  we  treat  a  uniform  half space  excited  by  a  one 

square  kilometer  fault,  centered  at  a  depth  of  five  kilometers.  The 

fault  is  vertical  strike  slip  and  the  source  time-function  is  shown 

in  Figure  1.  The  P  and  S  wave  speeds  (a  and  s)  and  density  (p)  of 

3 

the  half  space  are  6  km/sec,  3.46  km/sec,  and  2.7  gm/cm  , 
respectively.  Free  surface  displacements  computed  by  I4TRES  can  be 
compared  to  both  Cagniard-deHoop  and  finite  element  solutions.  The 
I4TRES  input  is  given  in  Table  B.l,  and  the  grid  configuration  is 
sketched  in  Figure  B.l. 

2.  TEST  PROBLEM  II:  BURIED  DISLOCATION  IN  A  LAYERED  ELASTIC 
HALF SPACE 

The  source  in  this  problem  is  identical  to  that  in  Problem  I. 
In  this  case,  however,  the  uniform  half space  is  overlain  by  two 

homogeneous  layers,  each  of  two  kilometers  thickness.  The  I4TRES 
input  is  given  in  Table  B.2,  the  material  properties  are  given  in 
Table  B.3,  and  the  grid  geometry  is  sketched  in  Figure  B.2. 

3.  TEST  PROBLEM  III:  SPONTANEOUS  RUPTURE  ON  A  PLANE 

In  this  problem,  rupture  is  forced  to  nucleate  in  a  circular 
region  of  a  fault  plane.  Subsequent  rupture  is  controlled  by  a 

slip-weakening  failure  criterion  (Ida,  1972;  Andrews,  1976).  In 

this  model,  a  finite  frictional  strength  aQ  is  assigned  to  the 

fault  plane  prior  to  initiation  of  sliding.  Slip  commences  when 
necessary  to  prevent  a  stress  concentration  in  excess  of  aQ  from 
occurring.  This  relative  displacement  is  assumed  to  weaken  the 
fault  plane,  as  in  Figure  B.3,  until  the  total  slip  equals  dQ,  at 


iL 


* 
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TABLE  B.1 

INPUT  FOR  I4TRES  TEST  PROBLEM  I 


♦ 


i 

! 


Grid  Limits 

JMAXG  =  60 
KMAXG  =  26 
LMAXG  =  51 

JPLMN 

JPLMX 

KPLMN  \  nQ  p]astiC  20nes 
KPLMX  1 
LPLMN  \ 

LPLMX  / 

JSLMN  =  6 

JSLMX  =  6 
KSLMN  =  1 
KSLMX  =  1 
LSLIP  =  26 


Zone  Sizes 

Ax  =  Ay  =  Az  =  1000  in  the  region  of  "uniform  zoning",  i.e., 
j  =  1  -  60,  K  =  1  -16,  L  =  11  -  41.  Outside  this  region,  the 
zone  size  increases  at  15  percent  per  zone  moving  away  from  the 
uniform  region. 


Total  time  steps  =  100 
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TABLE  B.l  (CONTINUED) 


Material  Properties 

Damping  factor 

BETA 

=  0.2 

Rotational  Q  factor 

BETAR 

=  1000 

Yield  stress 

YIELD 

=  1012 

Bulk  modulus 

AK 

=  5.4  x  101' 

Shear  modulus 

AM 

=  3.24  x  10 

Density 

RHO 

=  2700 

Source 

Prestress 

STLD 

=  108 

Sliding  friction 

SKIN 

=  0.9  x  108 

Rupture  velocity 

RUPV 

=  10® 

Hypocentral  node 

JFOCUS 

KFOCUS 

=  6 

=  1 

Boundary  Conditions 

J  =  1 

u. 

v,  w  all  free 

K  =  1 

u, 

w  fixed,  v  free 

L  =  1,  L  =  51,  J  =  60, 

K  =  26  u. 

v,  w  all  fixed 
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TABLE  B.2 


INPUT  FOR  I4TRES  TEST  PROBLEM  II. 


Grid  Limits 

JMAXG  =  60 
KMAXG  =  30 
LMAXG  =  59 


JPLMN 

JPLMX 

KPLMN 

KPLMX 

LPLMN 

LPLMX 


no  plastic  zones 


JSLMN  =  6 
JSLMX  =  6 
KSLMN  =  1 
KSLMX  =  1 
LSLIP  =  30 


Zone  Sizes 

Ax  =  Ay  =  Az  =  1000  in  the  region  of  "uniform  zoning,"  i.e., 
j  =  1  -  60,  k  =  1  -  16,  L  -  15  -  45.  Outside  this  region,  the 
zone  size  increases  at  15  percent  per  zone  moving  away  from  the 
uniform  region. 


Total  time  steps  =  150 
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TABLE  B.2  (CONTINUED) 


Material s 

Three  different  materials.  For  each  material,  the  damping  factors 
and  YIELD  will  be  BETA  =  0.2,  BETAR  »  1000,  YIELD  =  1012.  AK,  AMU 
and  ARHO  are  given  in  Table  B.3. 


Source 

STLD  =  108 
SKIN  =  0.9  x  108 
RUPV  =  108 
J FOCUS  =  6 
KFOCUS  =  1 


u,  v,  w  all  free 

u,  w  fixed,  v  free 

30  u,  v,  w  all  fixed 


Boundary  Conditions 
On  J  =  1 
On  K  =  1 

On  L  =  1,  L  =  59,  J  =  60,  K  = 


Output 

Saved  variables  are  u,  v,  w,  UD,  VD,  WD,  VSL,  VDSL. 
Nodes  at  which  output  is  required: 

(1,  1,  Z),  l  =  30  to  45 
(1 ,  k,  30),  k  =  1  to  16 
(1,  k,  29  +  k),  k  =  1  to  16 
(1,  5,  40) 

(1,  7,  44) 

!j  =  -1  to  1 

k  =  1  to  2 

Z  =  -1  to  1 
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30) 


30  plane  for  Test  Problem  II. 


SHEAR  STRESS 


SLIP 


Figure  B.3.  Slip-weakening  failure  model  (from  Andrews,  1976) 


which  state  the  fault  has  lost  all  cohesion  and  reaches  the  kinetic 
friction  level,  am¬ 
iable  B.4  gives  the  I4TRES  input  for  this  problem.  The  medium 
is  uniform  and  there  is  no  free  surface.  Test  problem  results  can 
be  compared  to  the  numerical  solution  obtained  with  the 

single-material  version  of  I4TRES.  These  results  are  reported  by 

Day  (1979,  Figures  3.13  and  3.14). 

4.  TEST  PROBLEM  IV:  BURIED  THRUST  FAULT 

In  this  test  case,  rupture  is  confined  to  a  small,  buried 
surface  to  which  a  free  surface  is  inclined  at  45  degrees.  The 
medium  is  a  uniform,  elastic  halfspace.  The  I4tres  input  is  given 
in  Table  B.5.  Figure  B.4  shows  the  configuration  of  the  inclined 
free  surface. 
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TABLE  B.4 

INPUT  FOR  I4TRES  TEST  PROBLEM  III 


Grid  Limits 

JMAXG  =  20 
KMAXG  =20 
LMAXG  =  41 

JPLMN 
JPLMX 
KPLMN 
KPLMX 
LPLMN 
LPLMX 

JSLMN  =  1 
JSLMX  =  10 
KSLMN  *  1 
KSLMX  =  10 
KSLIP  =  21 

Zone  Sizes 

Ax  =  Ay  =  Az  =  150  in  the  region  of  uniform  zoning,  i.e., 
j  =  1  -10,  k  =  1  -10,  £=12  -30.  Outside  this  region, 
the  zone  size  increases  at  15  percent  per  zone. 

Time  Stepping 


no  plastic  zones 


At  =  0.0125 

Total  time  steps  =  50 


TABLE  B.4  (CONTINUED) 


Material  Properties 


BETA 

= 

0.4 

BETAR 

= 

1500 

YIELD 

= 

1012 

AK 

= 

5.4  x  I0la 

AMU 

= 

3.24  x  1010 

RHO 

= 

2700 

Source 

Prestress 
Sliding  friction 
Strength  (a  ) 

Rupture  velocity 
Nucleation  radius 
1/2  dQ 

Hypocentral  node 

Boundary  Conditions 
J  =  1 
K  »  1 

L  *  1 ,  L  =  41 ,  J  =  20 


STLO 

= 

108 

SM1N 

0.9  x  108 

SMAX 

= 

1.02  x  108 

RUPV 

1732 

RCRIT 

= 

525 

SFRAC 

= 

.05 

J FOCUS 

= 

1 

KFOCUS 

= 

1 

u  fixed;  v. 

w  free 

u,  w  fixed; 

v  free 

II 

O 

u,  v,  w  all 

fixed 
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TABLE  B . 5 

INPUT  FOR  I4TRES  TEST  PROBLEM  IV 


Grid  Limits 


JMAXG  =  60 
KMAXG  =  39 
LMAXG  =  39 


JPLMN 

JPLMX 

KPLMN 

KPLMX 

LPLMN 

LPLMX 


no  plastic  zones 


JSLMN  =  I 
JSLMX  =  1 
KSLMN  =  17 
KSLMX  =  17 
LSLIP  =  17 

Zone  Sizes 

Ax  =  Ay  =  Az  =  1000  in  the  region  of  uniform  zoning,  i.e., 
3=1-  60,  k  =  11  -  29,  £  =  11  -  29.  Outside  this  region, 
the  zone  size  increases  at  15  percent  per  zone. 


At  =  0.08 

Total  time  steps  =  100 


TABLE  B. 5  (CONTINUED) 
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A  cross  section  of  the  grid  for  Test  Problem  IV,  showi 
the  free  surface  configuration. 
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APPENDIX  C 

CONVERSION  OF  TRES  TO  THE  CRAY  I,  AND  PARTIAL  VECTORIZATION 

In  this  appendix  we  record  our  experience  in  converting  the 
TRES  code  from  the  UNIVAC  1100/81  computer  to  the  CRAY  I  computer. 
This  discussion  is  included  in  hopes  of  assisting  other  prospective 
CRAY  users  in  estimating  the  effort  involved  in  code  conversion.  We 
emphasize  that  only  a  partial  vectorization  was  attempted  in  view  of 
the  available  resources  (i.e.,  2.5  man-months). 

We  discuss  first  the  initial  work  involved  in  getting  a 
working  code  without  vectorization.  The  original  TRES  code  was 
written  in  UNIVAC  FORTRAN  V,  which  includes  several  nonstandard 
statements  such  as  INCLUDE,  DEFINE,  RETURN  0,  and  some  special  I/O 
statements.  The  most  time-consuming  parts  of  the  initial  conversion 
process  were  conversion  of  input/output  operations  and  replacement 
of  INCLUDE  statements. 

We  will  start  wtih  the  input/output  operations  which  accounted 
for  most  of  the  difficulties  in  the  conversion  process.  The 
standard  FORTRAN  read  and  write  statements,  although  convenient  to 
use,  do  not  make  efficient  use  of  tapes  and  disks  because  the 
FORTRAN  language  does  not  permit  parallel  processing.  Furthermore, 
a  considerable  amount  of  time  is  used  in  processing  an  I/O  list 
because  of  its  generality.  To  overcome  those  problems,  the  UNIVAC 
came  out  with  nonstandard  FORTRAN  statements  wh’ch  allow  buffering. 
The  UNIVAC  1100/81  uses  NTRAN  I/O  package  which  provides  a  tool  for 
reading  and/or  writing  binary  information  on  tape  or  disk  files.  It 
also  provides  I/O  buffering  through  a  call  statement  in  the  FORTRAN 
language:  Call  NTRAN  (UNIT,  SEQUENCE  OF  OPERATION).  NTRAN  allows 
one  to  use  23  different  types  of  operations,  including  positioning 
on  the  disk  file.  The  hardware  architecture  of  the  UNIVAC  computer 
put  some  restrictions  on  the  FORTRAN  I/O  operations.  It  always 
reads  and  writes  a  multiple  of  28  words.  This  "MAGIC"  number  is 
used  to  calculate  the  new  position  on  the  disk.  We  also  have  found 


that  the  following  somewhat  unusual  situation  can  occur  and  is 
allowed  in  NTRAN:  one  can  write  initially  N  records  (N>1)  on  the 
disk,  rewind  the  file,  write  K  new  records  (K<N)  and  continue  to 
read  remaining  records. 

The  CRAY-1  FORTRAN  ( CFT )  also  allows  the  programmer  to  use 
buffering  operations,  but  its  repertoire  is  very  limited. 
Essentially,  only  five  operations  are  allowed:  BUFFER  IN,  BUFFER 
OUT,  BACK  SPACE,  BACKFILE  AND  ENDFILE.  To  perform  FORWARD  SPACING, 
one  has  to  use  BUFFER  IN  instruction.  Also,  the  READ  AFTER  WRITE 
combination  is  not  allowed  in  CFT. 

A  second  cumbersome  problem  was  the  replacement  of  INCLUDE 
statement,  which  is  used  to  insert  any  externally  defined  set  of 
UNIVAC  FORTRAN  statements  into  the  program  being  compiled.  The 
general  form  of  the  INCLUDE  statement  is  INCLUDE,  n,  where  n  is  the 
name  of  a  procedure  created  by  UNIVAC  procedure  definition  processor 
(PDP ) .  We  had  to  replace  all  INCLUDE  statements  in  each  subroutine 
by  set  of  PARAMETER,  TYPE,  and  DIMENSION  statements,  etc.  (Almost 
500  cards  had  to  be  punched  and  inserted  into  the  program  file). 
While  doing  these  replacements,  we  discovered  that  PARAMETER 
statement  in  the  two  compilers  are  slightly  different.  For  example, 
PARAMETER  N  =  3  on  the  UNIVAC  will  have  the  form  of  PARAMETER  (N  = 
3)  on  CRAY-I.  Also,  we  had  to  change  the  DEFINE  statements  used  in 
the  PDP  elements. 

One  also  can  expect  some  problems  working  with  Hollerith  and 
Alphanumeric  data  because  of  differences  in  the  word  length  (64  bit 
word  on  CRAY  versus  36  bit  word  on  UNIVAC  1100/81) 

Approximately  two  weeks  of  intensive  effort  by  a  programmer 
familiar  with  the  CRAY-1  were  required  to  obtain  a  diagnostic-free 
listing  of  the  CRAY  TRES  code.  Several  more  days  were  required  to 
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identify  and  correct  all  potential  divisions  by  zero,  a  circumstance 
tolerated  by  UNIVAC  (result  equals  zero)  and  exploited  in  the 
original  TRES  code  but  not  tolerated  by  the  CRAY.  Thus,  the  initial 
FORTRAN  conversion  phase  required  about  three  man-weeks. 

Next  we  discuss  the  vectorization  of  the  TRES  program.  We 
stress  that  our  goal  was  to  convert  the  existing  code  and  not  to 
fully  optimize  TRES,  which  would  require  rewriting  most  of  the  code 
and  is  considerally  beyond  the  scope  of  this  task.  The  first 
element  in  this  process  is  to  analyze  the  program  to  determine  where 
it  spends  its  time  and  to  concentrate  on  those  time-consuming 
parts.  We  have  found  that  in  our  case  almost  all  time  was  spent 
inside  SUBROUTINE  MOTION  and  some  subprograms  associated  with  it 
(ROTQ,  PLSTRN,  STRAIN).  The  original  code  is  overly  modular;  most 
of  the  subroutines  are  inside  the  DO  loops.  We  had  to  change  this 
structure  to  put  the  DO  loops  inside  the  subroutines.  This  done, 
our  efforts  were  directed  to  vectorization  of  MOTION,  PLSTRN,  etc. 
The  partial  vectorization  was  achieved  by  replacement  of  IF  loops  by 
DO  loops,  breaking  of  complicated  loops  to  several  short  and  simple 
ones,  removal  of  the  IF,  CALL  and  some  other  similar  statements  from 
the  DO-loops,  changing  the  order  of  the  nest  loops,  and  so  forth. 
The  CFT  diagnostic  is  very  helpful  in  this  matter,  providing  much  of 
the  needed  information  about  scalar  and  vector  loops. 

This  limited  vectorization  phase  of  the  conversion  task 
required  an  additional  three  man-weeks  for  completion.  The 
remainder  of  the  effort  on  the  task  was  expended  in  testing  and 
debugging. 


145 


